# How to monitor the progress of LinearSolve?

I am trying to do something like

d = RandomReal[{0, 50}, 3];
U = ResourceFunction["RandomMatrix"]["Orthogonal", Real, 3];
m = Transpose[U] . DiagonalMatrix[d] . U;
b = {6, 4, 5};
LinearSolve[N@m, N@b,
Method -> {"Krylov", Method -> "ConjugateGradient"}]


and want to monitor the decline of residue in the solving process. It seems that Solve and LinearSolve do not have the StepMonitor and EvaluationMonitor options, which are only available for iterative built-in functions.

I'm guessing blindly that, LinearSolve is not eligible for these options because even when an iterative method is specified, the iterating procedure is handled by low-level linalg libraries, rather than within the Wolfram language abstraction. If this is the case, it seems also to rule out the possibility of achieving this with Sow and Reap.

Is there some other way I have not think of that may do the trick?

Related posts:

Indirect workaround: Since this is not a real world problem, one can always use CG FindMinimum on the corresponding variational minimization problem, and StepMonitor the process. I however still appreciate further insights on the internals of LinearSolve. From the second linked post I understand that there is internally the function SparseArrayKrylovLinearSolve; is it possible to probe its implementation?

• Comment not directly related to your question: Instead of the resource function, one can use U = RandomVariate[CircularRealMatrixDistribution[3]]. Nov 15, 2022 at 9:01
• @user293787 Thanks! I intend to investigate the sampling measure of the resource function too Nov 15, 2022 at 9:38
• I'm guessing blindly that, LinearSolve is not eligible for these options That is correct. There is no direct way to monitor exact linear algebra computations. reference Daniel Lichtblau in comment to the first link. Nov 15, 2022 at 16:11
• I seem to have doomed myself with ever-prolonging time working on this, got to stop and take a note. How does built-in functions handle sub-options exactly? LinearSolve[m, b, Method -> {"Krylov", <suboptions>}] should take string option keys, but it also accepts symbol keys; is this intended? How can I do this in my custom functions? Nov 16, 2022 at 14:04

I think I got it. MSE rules and this should be obvious had I been competent.
This is what @unlikely did in his question

Reap@LinearSolve[m, b, Method -> {"Krylov",
"ResidualNormFunction" -> ((Sow[#, "residue_vector"]; Sow[Norm[#, 2]]) &)
}
]


This should intuitively suffice the purpose in question, however does not.

I finished my mock-up that I think is mathematically equivalent to @user293787 's

ClearAll[GravifermyKrylovLinearSolve]
Module[{symbol = GravifermyKrylovLinearSolve, $$symbol = SparseArrayKrylovLinearSolve}, ((MessageName[Evaluate@symbol, #] = MessageName[Evaluate@$$symbol, #]) &) /@ ToExpression[
FindList[
FileNameJoin[{SystemPrivate$$MessagesDir,$$Language,
"Messages.m"}], ToString[\$symbol]]
, StandardForm, HoldForm][[All, 1, 1, -1]]
];
Options[GravifermyKrylovLinearSolve] ^=
Options[SparseArrayKrylovLinearSolve];
SyntaxInformation[
GravifermyKrylovLinearSolve] ^= {"ArgumentsPattern" -> {_, _,
OptionsPattern[]},
"OptionNames" ->
Union[First /@ Options[GravifermyKrylovLinearSolve],
First /@ Options[SparseArrayKrylovLinearSolve]]};
GravifermyKrylovLinearSolve[m_, b_,
opt : OptionsPattern[SparseArrayKrylovLinearSolve]] :=
Module[{dim = Dimensions[m][[1]],
maxiter = OptionValue[MaxIterations], tol = OptionValue[Tolerance],
x0 = OptionValue["StartingVector"],
resnorm = OptionValue["ResidualNormFunction"] /. Automatic -> Norm,
CGfunc},

CGfunc =
Function[{n, A, \[FormalB], \[FormalCapitalN], \[Epsilon], x0},
Module[{x = x0, r = \[FormalB] - A . x0,
p = \[FormalB] - A . x0, \[Alpha], \[Beta]},
Do[If[resnorm[r] < \[Epsilon], Break[]];
If[p . A . p == 0,
Message[GravifermyKrylovLinearSolve::krystg]; Break[]];
\[Alpha] = Norm[r]^2/p . A . p;
x += \[Alpha] p; \[Beta] = Norm[r - A . (\[Alpha] p)]^2/
Norm[r]^2; r -= A . (\[Alpha] p); p = \[Beta] p + r,
\[FormalCapitalN]];
If[resnorm[r] > \[Epsilon],
Message[GravifermyKrylovLinearSolve::krymit, \
\[FormalCapitalN]]]; x
]];
If[Precision[{m, b}] == \[Infinity],
Message[GravifermyKrylovLinearSolve::kryinfp]];
(*
GravifermyKrylovLinearSolve::kryinz,"m"]];
GravifermyKrylovLinearSolve::kryinz,"b"]];
*)
Switch[OptionValue[Method],
CGfunc[dim, m, b, maxiter /. Automatic -> n,
tol /. Automatic -> 10^-12,
x0 /. Automatic -> ConstantArray[0, dim]]
]
] /; And[CheckArguments[GravifermyKrylovLinearSolve[m, b, opt], 2],
(If[MatrixQ[m], True,
Message[GravifermyKrylovLinearSolve::matrix, m, 1];
If[MatrixQ[Evaluate@m], True,
Message[GravifermyKrylovLinearSolve::krynfa, m, 1]]; False]),
(VectorQ[b] \[Or]
If[MatrixQ[b], True,
Message[GravifermyKrylovLinearSolve::matrix, b, 2]; False]),
((OptionValue[MaxIterations] == Automatic) \[Or]
If[IntegerQ@OptionValue[MaxIterations], True,
Message[GravifermyKrylovLinearSolve::kryit,
OptionValue[MaxIterations]]; False]),
((OptionValue["StartingVector"] == Automatic) \[Or]
If[VectorQ[OptionValue["StartingVector"], NumberQ], True,
Message[GravifermyKrylovLinearSolve::kryivec]; False]),
Switch[OptionValue[Method],
"SteepestDescent" | "Lanczos" | "Arnoldi" | "MinRES" | "GMRES" |
"CGR" | "CGNR" | "CGLE" | "QMR" | "BiCG" | "BiCGStab",
Message[GravifermyKrylovLinearSolve::krynimp]; False,
_, Message[GravifermyKrylovLinearSolve::kryme,
OptionValue["Method"]]; False
]
]


On the cases from his answer,

SeedRandom[1];
dim = 2000;
U = RandomVariate[CircularRealMatrixDistribution[dim]];
A = Transpose[U] . DiagonalMatrix[RandomReal[{1/100, 1}, dim]] . U;
b = RandomReal[{-1, 1}, dim];
{sol, reap} =
LinearSolve[A, b,
Method -> {"Krylov", Method -> "ConjugateGradient",
"ResidualNormFunction" -> ((Sow[Norm[#]]) &)}] // Reap;
{mysol, myreap} =
GravifermyKrylovLinearSolve[A, b, Method -> "ConjugateGradient",
"ResidualNormFunction" -> ((Sow[Norm[#]]) &)] //
Reap;
ListLinePlot[{reap[[1]], myreap[[1]]},
ScalingFunctions -> "Log"]
{sol, reap} =
LinearSolve[N@A, N@b,
Method -> {"Krylov", Method -> "ConjugateGradient",
"ResidualNormFunction" -> ((Sow[Norm[#]]) &)}] // Reap;
{mysol, myreap} =
"ResidualNormFunction" -> ((Sow[Norm[#]]) &)] //
Reap;
ListLinePlot[{reap[[1]], myreap[[1]]},
ScalingFunctions -> "Log"]


The results are

So maybe as @unlikely indeed suggested in his post, OptionValue["ResidualNormFunction"] is not called once per iteration for the built-in function.

• Good idea, but did you try it? If I try this using A and b in my answer, then the "ResidualNormFunction" seems to be evaluated only three times, and the only residue norms that one obtains are {44.7214,26.0726,1.43976*10^-8}. I do not think it only uses three steps? I use V12.3. Nov 16, 2022 at 4:56
• @user293787 I did some small scale tests, and I really assumed tha was the case. I'm a little confused right now... darn it Nov 16, 2022 at 7:08

Not directly an answer, but here is an implementation of conjugate gradient for a real, symmetric and positive definite matrix $$A$$. It Sows the norm of the residual at each step. This implementation uses more matrix-vector multiplications than necessary:

linearSolveConjugateGradient[A_,b_,tol_:10^(-10)] := Module[{p,x,r},
p = {};
x = 0.*b;
Do[
r = A.x-b;
If[Sow[Norm[r]]<tol,Break[]];
p = {r,First[p,Nothing]} // DeveloperToPackedArray;
p[[1]] = LeastSquares[p.A.Transpose[p],-p.r].p;
x += p[[1]];
,{Infinity}];
x
];


Example. Create a large (whatever that means) matrix $$A$$ with condition number about 100:

SeedRandom[1];
dim = 2000;
U = RandomVariate[CircularRealMatrixDistribution[dim]];
A = Transpose[U].DiagonalMatrix[RandomReal[{1/100,1},dim]].U;
b = RandomReal[{-1,1},dim];


Monitor progress:

ListLogPlot[Last[Reap[linearSolveConjugateGradient[A,b]]],
AxesLabel->{"step","Norm[residual]"}]


In this particular example, it is faster than the built-in conjugate gradient solver:

Norm[A.linearSolveConjugateGradient[A,b]-b] // AbsoluteTiming
(* {0.606085, 8.21629*10^-11} *)

`