# Algorithm for the resolvent of a matrix, and LUBackSubstitution

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?

• It might be fast to create a LinearSolveFunction using LinearSolve[T-z*IdentityMatrix[n]] since T is already upper triangular. May 27, 2022 at 17:41
• 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. May 27, 2022 at 18:25
• 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]}, LinearAlgebraBLASTRSV["U", "N", "N", tr, bb]; qq . bb] May 27, 2022 at 20:57
• 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. May 28, 2022 at 6:01
• Version 13.1 now has UpperTriangularMatrix[], which you can try for your resolvent computation. Jul 9, 2022 at 15:10

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];

• Warning: I have not been able to try any of this, since I am not on V13.1 yet. Jul 10, 2022 at 5:01
• user293787, you can always use the wolfram cloud to try out such things. Jul 10, 2022 at 23:55