2
$\begingroup$

I am computing U.U......U.x, where x is a vector, U a square matrix and the product involves N times U, i.e., (U^N).x if you like. One way to do this is

i=1;
While[i<=N,
x=U.x;
i=i+1;
]

Another way to do this, which saves time if N>>1 is really large, is to first multiply a few U matrices, say 10, and then let the loop run only N/10 many times:

U10=U.U.U.U.U.U.U.U.U.U;
i=1;
While[i<=N/10,
x=U10.x;
i=i+1;
]

Now, I made the observation that for large numerical U and x, the second method is significantly less accurate than the first. I assume the reason is that generating U10 introduces a non-negligible numerical error. Therefore, the question is:

How can I tell Mathematica to make the computation of U10 more accurate?

Thanks for your help!

$\endgroup$
2
  • 6
    $\begingroup$ Have you seen MatrixPower? $\endgroup$
    – Bob Hanlon
    Jan 19, 2023 at 0:20
  • $\begingroup$ Note N is a built-in, protected function. $\endgroup$
    – Michael E2
    Jan 19, 2023 at 1:14

1 Answer 1

3
$\begingroup$
Clear["Global`*"]

SeedRandom[1234];

u = RandomReal[10, {3, 3}];

To increase the precision of the calculations, use arbitrary precision rather than machine precision.

u = u // SetPrecision[#, 15] &;

n = 10;

t1 = RepeatedTiming[sol1 = MatrixPower[u, n]][[1]]

(* 0.0000353723 *)

t2 = RepeatedTiming[sol2 = Nest[# . u &, u, n - 1]][[1]]

(* 0.00011738 *)

MatrixPower is much more efficient

t2/t1

(* 3.31842 *)

The results are Equal

sol1 == sol2

(* True *)

There was only a limited loss of precision from the calculations using arbitrary precision

Precision /@ {sol1, sol2}

(* {14., 14.} *)
$\endgroup$
4
  • 5
    $\begingroup$ One reason you have minimal precision loss is that all the entries in u are positive (no subtractive cancellation). Using u = RandomReal[{-2, 2}, {3, 3}] and n = 50, I get a rather significant difference between MatrixPower and Nest (MatrixPower being the better of the two). $\endgroup$
    – Michael E2
    Jan 19, 2023 at 2:05
  • $\begingroup$ Thanks for pointing this out, but it does not really answer my question (though it speeds up the code a bit, which is helpful). However, the problem is still that U10.x with U10 = MatrixPower[U,10] gives a DIFFERENT result from applying U iteratively 10 times directly to x. It is precisely this issue that I like to get settled (if possible). Increasing precision seems not to be an option as it slows down the code a lot... $\endgroup$ Jan 19, 2023 at 23:06
  • 1
    $\begingroup$ If you are having precision issues then you must track and control the precision, i.e., use arbitrary precision rather than machine precision. In my example, both methods gave the same result to within allowed tolerances. You should provide a concrete example that demonstrates the issue that you are having. $\endgroup$
    – Bob Hanlon
    Jan 19, 2023 at 23:25
  • $\begingroup$ Let's say U = U[dt] is some approximate "time-evolution operator" depending on a step size dt. Now, the numerics were converged for some dt by applying U[dt] ten times to some vector x. However, first computing U[dt]^10 and applying this to the same vector showed convergence only for a smaller dt! I think the only explanation is that Mathematica introduces a smaller numerical error when multiplying a matrix with a vector than when multiplying a matrix with a matrix (even if the final result is still applied to the same vector). Thanks for helping me to sort this out! $\endgroup$ Jan 20, 2023 at 0:14

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