3
$\begingroup$

I have $N\ 2 \times 2$ matrices, each one containing one variable energy which may be complex. I want to muliply them all together to a single matrix M and then find all the energy for which M[[2,1]] + M[[2,2]] = 0.

What I do so far is the following:

I generate the matrices and put them all inside a list called matrixlist. Each one still contains the variable 'energy'.

M = Fold[Dot, IdentityMatrix[2], matrixlist]; 
ev = FindRoot[M[[2, 1]] + M[[2, 2]] == 0, {energy, 1.1}]] 

This prodecure sometimes finds one of the complex energy.

The problem I'm having is with memory. As soon as I multiply more than, say, 20 matrices, my PC gives up (kernel shutdown). I would like to handle $N= 100$.

Any idea how to find all the complex solutions energy within reasonable time? Maybe it helps for you to know that all the matrices are unitary (but all are different and the energy dependency is complicated)?

edit: sorry for the mistakes, I am new to SE, thank you for welcoming me! Here is the relevant code:

   V = x^2 - 0.1 x^4 I - 0.1 I + 1;  
h = 1;
dischalf = Table[{x - h/2, V}, {x, 0, 10, h}];
k[j_] := Sqrt[2 (energy - U[j])];
U[j_] := dischalf[[j + 1, 2]];
B[j_, i_] := 
  1/(2 k[j]) {{(k[j] + k[i]) Exp[
       I (k[j] + k[i]) h/2], (k[j] - k[i]) Exp[
       I (k[j] - k[i]) h/2]}, {(k[j] - 
        k[i]) Exp[-I (k[j] - k[i]) h/2], (k[j] + 
        k[i]) Exp[-I (k[j] + k[i]) h/2]}}; 
SmatrixlistHalb = 
  Reverse[Table[B[j + 1, j], {j, 0, Length[dischalf] - 2}]];
STransferMatrixH = Fold[Dot, IdentityMatrix[2], SmatrixlistHalb];
Timing[FindRoot[
  STransferMatrixH[[2, 1]] + STransferMatrixH[[2, 2]] == 0, {energy, 
   0.5}]]

This calculation takes on my PC ~15sec. But I would like to decrease the stepsize h to 0.1, which is a bit too much for my good old calculator.
To see better how the matrices look like, here the matrix $B[1,2]$: $B[1,2]= \left( \begin{array}{cc} \frac{e^{\frac{1}{2} i h (k(1)+k(2))} (k(1)+k(2))}{2 k(1)} & \frac{e^{\frac{1}{2} i h (k(1)-k(2))} (k(1)-k(2))}{2 k(1)} \\ \frac{e^{-\frac{1}{2} i h (k(1)-k(2))} (k(1)-k(2))}{2 k(1)} & \frac{e^{-\frac{1}{2} i h (k(1)+k(2))} (k(1)+k(2))}{2 k(1)} \\ \end{array} \right)$

If I can help you helping me in any way let me know!

$\endgroup$
  • 1
    $\begingroup$ Would you edit into your post a half dozen representative matricies from your matrixlist so that readers might have an idea what those look like? Does M=Dot@@matrixlist give you exactly the same result as your M=Fold[...] does (for a half dozen matricies in matrixlist)? Does your kernel fail during the Fold or during the FindRoot? $\endgroup$ – Bill Jan 13 '16 at 21:14
  • 2
    $\begingroup$ E is a reserved symbol denoting the base of the natural log./exp. functions. $\endgroup$ – Michael E2 Jan 13 '16 at 21:17
  • $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$ – Michael E2 Jan 13 '16 at 21:18
  • $\begingroup$ In general, avoid starting variables with capital letters to avoid clashing with built-in functions. $\endgroup$ – N.J.Evans Jan 13 '16 at 21:18
  • $\begingroup$ Maybe Expand[Dot[##]] & will do better than Dot alone. $\endgroup$ – Daniel Lichtblau Jan 13 '16 at 21:45
4
$\begingroup$

General idea

A general rule of thumb to reduce the time and memory needed for a computation is to use inexact (floating-point) numerics instead of exact or symbolic approaches, and to use it early. Of course, one might worry about accuracy, but there's no guarantee that feeding an exact, symbolic problem to FindRoot will be more accurate than feeding an inexact one.

For a problem $f(x) = 0$, from $dy = f'(x_0) \; dx$, the norm of the error $|dy|$ in the residual $f(x_0^*)$ of best floating-point approximation $x_0^*$ to an exact solution $x_0$ will roughly at most the norm of the Jacobian (derivative) $|f'(x_0)|$ times the uncertainty $|dx|$ in $x_0$ from the floating-point approximation. The uncertainty $|dx|$ is at most 0.5 ulp (unit in the last place) and is approximately $x_0\varepsilon/2$, where $\varepsilon$ is $MachineEpsilon. So a residual that satisfies $$f(x_0) \le |f'(x_0^*)|\, x_0\,\varepsilon/2$$ represents a solution that is probably pretty good.

By this measure, for the example problem, machine precision does a pretty good job, so compiling some of the computation will speed things up more. Otherwise, one could use arbitrary-precision numbers and replaced the compiled functions with regular Mathematica functions.

OP's code

The OP's code takes 10+ seconds and seems to have a large residual but not relative to the magnitude of the derivative (~ 10^51).

AbsoluteTiming[
 FindRoot[STransferMatrixH[[2, 1]] + STransferMatrixH[[2, 2]] == 
   0, {energy, 0.5}]]
STransferMatrixH[[2, 1]] + STransferMatrixH[[2, 2]] /. Last[%]

FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances. >>

(*
  {10.7866, {energy -> 1.78929 - 0.154278 I}}
  4.15384*10^35 + 9.96921*10^35 I
*)

Speed up: Compile the objective function

We compile the objective function for speed, objCF, and provide a _?NumericQ-protected interfaces, obj, so that the compiled functions won't evaluate symbolically. (I had done the same for the Jacobian, but it turns out not to help that much.)

Clear[obj, jac];
With[{mat = Dot @@ SmatrixlistHalb},  (* This takes a long time if  h  is small *)
  objCF = Compile[{{energy, _Complex}},
    Evaluate[mat[[2, 1]] + mat[[2, 2]]]
    ];
  obj[e0_?NumericQ] := objCF[e0];
  jacCF = Compile[{{energy, _Complex}},
    Evaluate@D[{mat[[2, 1]] + mat[[2, 2]]}, {{energy}}]];
  jac[e0_?NumericQ] := jacCF[e0];
  ];

The objective function alone give the same result, albeit in much less time:

(sol = FindRoot[obj[e0] == 0, {e0, 1/2}]) // AbsoluteTiming
obj[e0] /. sol

FindRoot::lstol: ...

(*
  {0.045566, {e0 -> 1.78929 - 0.154278 I}}
  -6.64614*10^35 + 6.64614*10^35 I
*)

By supplying the Jacobian, the root-finding is faster and more accurate:

(sol = FindRoot[obj[e0] == 0, {e0, 0.5}, Jacobian -> jac[e0]]) // AbsoluteTiming
obj[e0] /. sol

FindRoot::lstol: ...

(*
  {0.046047, {e0 -> 1.78929 - 0.154278 I}}
  -5.81537*10^35 + 3.32307*10^35 I
*)

Using less memory

Computing the matrices for a given value of energy and then their product saves on memory and it doesn't take any longer. Again I tried with the Jacobian, but again it didn't help. Code note: The form Fold[Dot, matCF[e0]] work in V10 and above; for earlier versions of Mathematica, use the OP's Fold[Dot, IdentityMatrix[2], matCF[e0]].

Clear[obj2, jac2];
With[{
   matlist = SmatrixlistHalb,
   dmatlist = D[SmatrixlistHalb, energy]},
  matCF = Compile[{{energy, _Complex}}, matlist];
  dmatCF = Compile[{{energy, _Complex}}, dmatlist];
  obj2[e0_?NumericQ] := With[{mat2 = Fold[Dot, matCF[e0]]},
    mat2[[2, 1]] + mat2[[2, 2]]
    ];
  jac2[e0_?NumericQ] := With[{m = matCF[e0], dm = dmatCF[e0]},
    With[{j = Sum[Fold[Dot, ReplacePart[m, i -> dm[[i]]]], {i, Length@m}]},
     {{j[[2, 1]] + j[[2, 2]]}}]];
  ];

We get similar results:

(sol = FindRoot[obj2[e0] == 0, {e0, 1/2}]) // AbsoluteTiming
obj2[e0] /. sol

FindRoot::lstol: ...

(*
  {0.047256, {e0 -> 1.78929 - 0.154278 I}}
  1.24615*10^35 + 0. I
*)

Numerical stability

The large magnitude of the Jacobian at the solution,

absJ = Abs@ jac[e0] /. {e0 -> 1.7892907120093438` - 0.15427772564884512` I}
absDX = Abs[e0] /. {e0 -> 1.7892907120093438` - 0.15427772564884512` I}
(* 
  {{4.04292*10^51}}
  1.79593
*)

makes it hard to get an exact root at machine precision. One should expect an a residual of approximately at most

absJ * absDX * $MachineEpsilon / 2
(*  {{8.06111*10^35}}  *)

which is what we see in the results of FindRoot. Thus it is a happy accident that obj has an apparently exact root, but obj2 does not:

obj[e0] /. {e0 -> 1.7892907120093438` - 0.15427772564884512` I}
obj2[e0] /. {e0 -> 1.7892907120093438` - 0.15427772564884512` I}
(*
  0. + 0. I
  2.4923*10^35 + 3.32307*10^35 I
*)
$\endgroup$
  • $\begingroup$ Thank you alot for giving such a detailed answer! I do not understand most of the code, so I just copy and pasted it. I will go through it line by line and use the documentary to fully understand what you are doing. Coying and pasting it I noticed two things: 1. My version of mathematica (v. 10) tells me something is missing in the 2 Fold commands in obj2[], inserting a IdentityMatrix[2] right after "Dot" fixes the problem. 2 The calcuation for h=1, is incredibe much faster, thank you alot for that. However, if I half the stepsize (h=0.5), both, obj[e0] and obj2 give me a kernel shutdow $\endgroup$ – Luke Jan 14 '16 at 13:36
  • $\begingroup$ @Luke You're welcome. 1. In V10, Fold[f, list] is equivalent to Fold[f, First[list], Rest[list]], but in V9 and earlier you have to write it explicitly as 3 arguments (or use the identity matrix as you do). (This was faster for me than Dot @@ matlist, when matlist was a "packed array," which you can search the site for an explanation.) 2. Use obj2 and omit mat = Dot @@ SmatrixlistHalb, which is unnecessary and I will delete presently. $\endgroup$ – Michael E2 Jan 14 '16 at 13:55
  • $\begingroup$ wow it works perfectly! Thank you, thank you! $\endgroup$ – Luke Jan 14 '16 at 14:40
  • $\begingroup$ @Luke You're welcome. Let me know if you have questions. I should point out that I didn't really end up using the derivative, except in the analysis, so dmatlist = D[SmatrixlistHalb, energy]}, dmatCF = Compile[{{energy, _Complex}}, dmatlist], and jac2[e0_?NumericQ] := ... could be omitted. $\endgroup$ – Michael E2 Jan 14 '16 at 17:33
  • $\begingroup$ I want to let you know, that I am back for more questions :) The new question is quite lengthy so i postet it here: mathematica.stackexchange.com/questions/107603/… . I was wondering if maybe your jacobi approach might help. I did not completely understand what it does so I have some troubles implementing it myself... $\endgroup$ – Luke Feb 17 '16 at 17:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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