2
$\begingroup$

Hello I'm trying to compute the trace of the product of matrices like this

$$\begin{equation} \begin{pmatrix} a+\imath b &c &0 \\ c& d &c \\ 0& c &f+\imath g \end{pmatrix}^{-1} \begin{pmatrix} -2b &0 &0 \\ 0 & 0 &0 \\ 0& 0 &0 \end{pmatrix} \begin{pmatrix} a+\imath b &c &0 \\ c& d &c \\ 0& c &f+\imath g \end{pmatrix}^{*-1} \begin{pmatrix} 0 &0 &0 \\ 0 & 0 &0 \\ 0& 0 &-2g \end{pmatrix} \end{equation}$$

where $*$ is the comlpex transpose and all the components are real.

This situation is very easy to compute it but if I want to do it for a 50x50 matrix (the first is a tridiagonal matrix with $d$ in the main diagonal except the $(1,1)$ and $(N,N)$ elements, $c$ in the other two diagonals,the second and the forth have non zero elements only in the $(1,1)$ and $(N,N)$ positions) it's almost take a lot of time.

Taking advantage of the symmetry and that at the end only the last column is non zero (and therefore only the element $(N,N)$ will contribute) is there a way to compute it faster?

Thanks

Edit: The code I'm using for a 10x10 matrix is the following

    M = {{a + I b, c, 0, 0, 0, 0, 0, 0, 0, 0}, {c, d, c, 0, 0, 0, 0, 0, 0,
 0}, {0, c, d, c, 0, 0, 0, 0, 0, 0}, {0, 0, c, d, c, 0, 0, 0, 0, 
0}, {0, 0, 0, c, d, c, 0, 0, 0, 0}, {0, 0, 0, 0, c, d, c, 0, 0, 
0}, {0, 0, 0, 0, 0, c, d, c, 0, 0}, {0, 0, 0, 0, 0, 0, c, d, c, 
0}, {0, 0, 0, 0, 0, 0, 0, c, d, c}, {0, 0, 0, 0, 0, 0, 0, 0, c, 
f + I g}}; 

    G1 = {{-2 b, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 
0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0,
 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 
0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0,
 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 
0}};
    G2 = {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 
0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0,
 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 
0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 
0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 
0, -2 g}};

    Inverse[M].G1.Inverse[ConjugateTranspose[M]].G2 //ComplexExpand

To make tridiagonal matrices I'm using the comand

    SparseArray[{Band[{1, 2}] -> c, Band[{1, 1}] -> d, Band[{2, 1}] -> c},10] // Normal

and then recopying it, change only the two elements I need to change and then I operate

Improve this question
$\endgroup$
7
  • $\begingroup$ Do you want a symbolic answer in terms of $a,b,c,d,f,g$, or are you happy with a function that takes numerical values for them and returns the trace? The two will obviously perform rather differently. $\endgroup$
    – Emilio Pisanty
    Commented Aug 12, 2016 at 14:31
  • $\begingroup$ @EmilioPisanty symbolic answer because each of this terms are functions of other variables but always this varables appear on the same way $\endgroup$
    – Daniel
    Commented Aug 12, 2016 at 14:35
  • $\begingroup$ Can you enter your matrices in copy-and-paste-able Mathematica code, properly formatted in code blocks? It will make it much easier for use to help you if we can just copy and paste into our own copies of Mathematica rather than have to type it in ourselves. $\endgroup$
    – march
    Commented Aug 12, 2016 at 15:35
  • $\begingroup$ @march I edited the post with the code I'm using for a 10x10 matrix $\endgroup$
    – Daniel
    Commented Aug 12, 2016 at 15:45
  • $\begingroup$ I'll think about this, but it might be worth writing this using the Hadamard product (i.e. element-wise multiplication of matrices) in this form. $\endgroup$
    – march
    Commented Aug 12, 2016 at 15:55

1 Answer 1

Reset to default
1
$\begingroup$

The following exploits sparse matrices, and uses LinearSolve to avoid computing the inverse which is not sparse. I suspect that further optimisation would be possible, but it is certainly practical for problems of the scale you describe. (I have assumed that you requre a numeric, rather than a symbolic solution)

test[n_] :=  Module[{c = 1.3, b = 0.2, a = 4.1, f = 2.2, g = -0.3, d = 0.5, M, G1, G2},
M = SparseArray[{Band[{1, 2}] -> c, Band[{1, 1}] -> d, Band[{2, 1}] -> c}, n];
M[[1, 1]] = a + I b;
M[[-1, -1]] = f + I g;
G1 = SparseArray[{1, 1} -> -2 b, {n, n}];
G2 = SparseArray[{-1, -1} -> -2 g, {n, n}];
Tr[LinearSolve[M, G1].LinearSolve[ConjugateTranspose[M], G2]]]

In my tests, this agrees with the formula that you gave. I have

AbsoluteTiming[test[5000]]
(* {9.70897, -0.0233494 + 7.80626*10^-17 I} *)
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.