Kronecker Delta functions and block matrices

[Andrei Berceanu]

Hi guys,

I am trying to programatically generate a Julia function that calculates a 6x6 matrix.
I have written a short LaTeX description of the problem at:

http://nbviewer.ipython.org/github/berceanu/notebooks/blob/master/OPODrag/block_matrix.ipynb

The basic idea is, I have an analytical expression for the matrix elements of the 3x3 matrices M and Q, M_{m,n} representing the element on row m, column n. The problem is that these depend on sums over Kronecker symbols, which need to be contracted. Once I build the M and Q matrices, I need to join them in the block-form given in the above notebook.
A working strategy that I have so far is to generate the matrix L in Mathematica, and copy-paste the resulting (huge) expression to Julia. I would now like to bypass Mathematica completely, and thought of using Sympy. While it has the capability of contracting the sums over the delta symbols, it is not clear to me how to then export the resulting function to Julia. If this can be done without Sympy, i.e. in pure Julia, that would be awesome!

On a last note, performance of the resulting function is important, because I will need to digonalize a lot of these 6x6 matrices, calculated for different parameter values (eg. varying k_p, k_s, omega_p, etc).

[Steven G. Johnson]

It should be straightforward to write a macro that simplifies the Kronecker sums to give you an efficient Julia expression.

[Tomas Lycken]

Maybe https://github.com/Jutho/TensorOperations.jl can be helpful here? I've never used it myself, but just by the name and the description ("Julia package for tensor contractions and related operations") it sounds like exactly what you're looking for...

// T

[Andrei Berceanu]

Steven, could you detail your macro proposal please? What exactly do you have in mind?

[Steven G. Johnson]

On Thursday, August 28, 2014 6:36:16 AM UTC-4, Andrei Berceanu wrote:

Steven, could you detail your macro proposal please? What exactly do you have in mind?

e.g. you could have a macro like (caution: untested):

macro M(m,n, ε1, ε2, ε3, X1,X2,X3, A1,A2,A3)
   ε = [ε1, ε2, ε3]
   X = [X1,X2,X3]
   A = [A1,A2,A3]
   Asum = {:+}
   for q=1:3, t=1:3
     if m+q == n+t
         push!(Asum, :($(A[q])' * $(A[t])))
     end
   end
   if m == n
      :($(ε[m]) + 2*abs2($(X[m])) * $(Expr(:call, Asum)))
   else
      :(2 * $(X[m])' $(X[n]) * $(Expr(:call, Asum)))
   end
end

Then calling e.g. @M(1,2,ε1, ε2, ε3, X1,X2,X3, A1,A2,A3) would insert a hard-coded
expression for M[1,2] according to your formula, with only the nonzero delta terms.  Here εm denotes the m'th case of your expression [ε(k_m+k) .... ] in brackets, Xm denotes X(k_m+k), and Am denotes Ã_m.

Macros are essentially functions evaluated at parse-time, not at runtime, so they are exactly what you want in order to auto-generate fast inlined code resulting from some complicated symbolic expression that is known in advance, with whatever simplifications you can devise.

You could write another macro for the entries of Q, and another macro to put it all together and evaluate the whole matrix if you want.

[Andrei Berceanu]

Thanks a lot for the detailed answer. I changed a bit your macro, introducing conj() instead of ' and splicing the Kronecker sum to get a functional version:

macro M(m,n, ε1,ε2,ε3, X1,X2,X3, A1,A2,A3)
   ε = [ε1,ε2,ε3]

   X = [X1,X2,X3]
   A = [A1,A2,A3]
   Asum = {:+}
   for q=1:3, t=1:3
     if m+q == n+t

         push!(Asum, :(conj($(A[q])) * $(A[t])))

     end
   end
   if m == n

      :($(ε[m]) + 2*abs2($(X[m])) * $(Expr(:call, Asum...)))
   else
      :(2 * conj($(X[m])) * $(X[n]) * $(Expr(:call, Asum...)))
   end
end

I also defined a macro for the elements of Q:
macro Q(m,n, X1,X2,X3, A1,A2,A3)

   X = [X1,X2,X3]
   A = [A1,A2,A3]
   Asum = {:+}
   for q=1:3, t=1:3

     if m+n == q+t

         push!(Asum, :($(A[q]) * $(A[t])))
     end
   end

    :(conj($(X[m])) * conj($(X[n])) * $(Expr(:call, Asum...)))
end

and now I'm a bit stuck on putting it all together to get the full L matrix. Im trying to calculated L(i,j) based on M and Q for the 4 separate blocks. I only implemented the upper left block but I get an error

macro L(i,j, ε1,ε2,ε3, X1,X2,X3, A1,A2,A3)
    if i<4
        if j<4
            @M(i,j, ε1,ε2,ε3, X1,X2,X3, A1,A2,A3)
        end
    end
end

`+` has no method matching +(::Symbol, ::Int64)

It seems that the @M macro needs a value for the m,n indices for the comparison
if m+q == n+t

Any ideas how to solve this? The end-result I want is a function L(k, A1, A2, A3, etc) which would return the full 6x6 matrix.

//A

[Steven G. Johnson]

I think you want

quote
     @M($i,$j, $ε1, ....)
end

in order for your L macro to return an expression (that in turn calls the M macro).

Macros are not functions.  When you call @M(i,j,...), it doesn't pass the *value* of the arguments i,j, etcetera, it passes the symbolic expressions :i, :j, and so on.

Of course, if you are only calling M from inside another macro, it is perfectly valid to just change M from a macro to an ordinary function (returning an expression) ... just change "macro" to "function" in the definition of M, and call it without the "@" inside L.

[Andrei Berceanu]

This is what I came up with for the full matrix calculation:

https://github.com/berceanu/notebooks/blob/master/julia/macros.ipynb

I now have a macro for the elements of M, another for the elements of Q and yet another for the elements of L. In the end I define the genmatL function which loops over the dimensions of L and calculates all the matrix elements. Does this look sound to you?

A few things that are still unclear:
1. Can this approach be optimized further? I can't help the feeling that it contains a fair amount of redundancy, especially looking at the quote $var end kind of expressions. I can replace the macros with functions as you suggest, but then it's not clear to me when to use a function and when a macro.
2. I would like to see the functional form of the full matrix L (not just element by element, but the whole matrix of expressions), in terms of the functions \epsilon, \gamma and X. I tried with macroexpand acting on genmatL, but that doesn't seem to work. This is basically just to convince myself that the approach indeed works as intended.

Thanks again!

[Andrei Berceanu]

On a more general note, I am thinking that with a bit of cleanup, this code could make a good Julia metaprogramming example, especially since I haven't seen many of those around :)

//A

[Steven G. Johnson]

On Saturday, August 30, 2014 4:52:49 AM UTC-4, Andrei Berceanu wrote:

A few things that are still unclear:
1. Can this approach be optimized further? I can't help the feeling that it contains a fair amount of redundancy, especially looking at the quote $var end kind of expressions. I can replace the macros with functions as you suggest, but then it's not clear to me when to use a function and when a macro.

There's no practical difference between a macro that calls a function (at macro-evaluation time, not in the generated expression) and a macro that calls a macro

What I usually do is to write *everything* as functions first, to make it easier to look at the resulting expressions.  Once the resulting expressions are what I want, then I convert the top-level function into a macro.

In terms of optimization, the speed of the macro is irrelevant; the question is the efficiency of the generated code.  Look at the output expression and see if you can identify opportunities for further optimization.
 

2. I would like to see the functional form of the full matrix L (not just element by element, but the whole matrix of expressions), in terms of the functions \epsilon, \gamma and X. I tried with macroexpand acting on genmatL, but that doesn't seem to work. This is basically just to convince myself that the approach indeed works as intended.

"That doesn't work" is hard to respond to.   Macros can do any code generation you want; it's just a matter of debugging.

Like I said, I would write things as functions first (and pass symbolic expressions for the arguments A1,A2,A3 etcetera) and tweak it until the resulting expression is what you want, before converting the top-level function to a macro.