6
$\begingroup$

Suppose A is a complex n x n matrix where n is quite big, given as a numeric array. For many values of the complex number z, and many complex vectors b with n components, I need to compute LinearSolve[A-z*IdentityMatrix[n],b]. But calling LinearSolve repeatedly is slow, and I need something faster. Note that the matrix Inverse[A-z*IdentityMatrix[n]] is also known as the resolvent of A.

In my application, A is not diagonalizable, so I cannot make an eigendecomposition. But I can always assume that z is not an eigenvalue of A. I think a reasonable algorithm is to use a Schur decomposition as a pre-processing step, and to perform a backsubstitution each time the resolvent is applied:

(* random example *)
n = 2000;
A = RandomComplex[{-1-I,1+I},{n,n}];

(* slow implementation *)
resolventSlow[z_,b_] := LinearSolve[A-z*IdentityMatrix[n],b];

(* preprocessing *)
{Q,T} = SchurDecomposition[A,RealBlockDiagonalForm->False];

(* faster implementation *)
resolventFast[z_,b_] := Q.LUBackSubstitution[{T-z*IdentityMatrix[n],Range[1,n],1.},Conjugate[Conjugate[b].Q]];

(* timing *)
z = RandomComplex[{-1-I,1+I}];
b = RandomComplex[{-1-I,1+I},n];
RepeatedTiming[x1 = resolventSlow[z,b];] (* about 0.5 seconds *)
RepeatedTiming[x2 = resolventFast[z,b];] (* about 0.1 seconds *)
Chop[Norm[x1-x2]] (* zero as expected *)

The function LUBackSubstitution that I use for backsubstitution has been obsolete since 2003, and its specification is a little arcane. See for example this v4 documentation and this documentation of LUDecomposition.

Question: Is there an idiomatic replacement for LUBackSubstitution in the code above, that does not explicitly invoke obsolete symbols or low-level libraries such as BLAS? Alternatively, is there another useful algorithm for the resolvent that I could use instead?

$\endgroup$
8
  • 2
    $\begingroup$ It might be fast to create a LinearSolveFunction using LinearSolve[T-z*IdentityMatrix[n]] since T is already upper triangular. $\endgroup$ May 27, 2022 at 17:41
  • 1
    $\begingroup$ Ouch. Yes, I had to see for myself after you noted this, and that is pretty bad. I do not have an alternative to propose though. That stated, I don't think LUBackSubstitution is going anywhere so it should be safe to use, deprecation notwithstanding. $\endgroup$ May 27, 2022 at 18:25
  • 4
    $\begingroup$ It's unfortunate that you give the restriction "does not explicitly invoke ... low-level libraries such as BLAS"; otherwise: {qq, tt} = SchurDecomposition[A, RealBlockDiagonalForm -> False]; resolventFast2[z_, b_] := Module[{tr = tt - z IdentityMatrix[n], bb = b . Conjugate[qq]}, LinearAlgebra`BLAS`TRSV["U", "N", "N", tr, bb]; qq . bb] $\endgroup$ May 27, 2022 at 20:57
  • 1
    $\begingroup$ Thank you very much, your TRSV code works for me and is about as fast as the LUBackSubstitution code. I explicitly excluded BLAS since I was hoping there could be, for example, some way to give LinearSolve a hint about upper triangularity. $\endgroup$
    – user293787
    May 28, 2022 at 6:01
  • 1
    $\begingroup$ Version 13.1 now has UpperTriangularMatrix[], which you can try for your resolvent computation. $\endgroup$ Jul 9, 2022 at 15:10

1 Answer 1

3
$\begingroup$

As pointed out by @J.M., starting with Version 13.1 (June 2022) there is UpperTriangularMatrix. It is efficient with LinearSolve, see the Details-section of the documentation.

Therefore one can now use the following idiomatic code:

(* V13.1 *)

(* preprocessing *)
{Q,T} = SchurDecomposition[A,RealBlockDiagonalForm->False];

(* new solution *)
resolventFast[z_,b_] := Q.LinearSolve[UpperTriangularMatrix[T-z*IdentityMatrix[n]],
                                      Conjugate[Conjugate[b].Q]];

Possible improvement 1. To also exploit efficient storage, one could call T = UpperTriangularMatrix[T] right after SchurDecomposition.

Possible improvement 2. Version 13.1 also has BlockUpperTriangularMatrix (experimental). If A is real then one can use SchurDecomposition in default mode which produces a block-upper-triangular T, but has the advantage that both Q and T are real which saves memory:

(* V13.1 -- only for real A *)

{Q,T} = SchurDecomposition[A];
resolventFast[z_,b_] := Q.LinearSolve[BlockUpperTriangularMatrix[
                                        T-z*IdentityMatrix[n]],b.Q];
$\endgroup$
2
  • 1
    $\begingroup$ Warning: I have not been able to try any of this, since I am not on V13.1 yet. $\endgroup$
    – user293787
    Jul 10, 2022 at 5:01
  • 1
    $\begingroup$ user293787, you can always use the wolfram cloud to try out such things. $\endgroup$ Jul 10, 2022 at 23:55

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.