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:

  1. Monitoring the execution of Solve - @StarBucK
  2. What does LinearSolve Preconditioner do? - Answer by @Henrik Schumacher

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 SparseArray`KrylovLinearSolve; is it possible to probe its implementation?

  • $\begingroup$ Comment not directly related to your question: Instead of the resource function, one can use U = RandomVariate[CircularRealMatrixDistribution[3]]. $\endgroup$
    – user293787
    Nov 15, 2022 at 9:01
  • $\begingroup$ @user293787 Thanks! I intend to investigate the sampling measure of the resource function too $\endgroup$
    – Gravifer
    Nov 15, 2022 at 9:38
  • 3
    $\begingroup$ 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. $\endgroup$
    – Nasser
    Nov 15, 2022 at 16:11
  • $\begingroup$ 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? $\endgroup$
    – Gravifer
    Nov 16, 2022 at 14:04

2 Answers 2


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

Module[{symbol = Gravifer`myKrylovLinearSolve, $symbol = 
  ((MessageName[Evaluate@symbol, #] = 
       MessageName[Evaluate@$symbol, #]) &) /@ ToExpression[
      FileNameJoin[{System`Private`$MessagesDir, $Language, 
        "Messages.m"}], ToString[$symbol]]
     , StandardForm, HoldForm][[All, 1, 1, -1]]
Options[Gravifer`myKrylovLinearSolve] ^= 
   Gravifer`myKrylovLinearSolve] ^= {"ArgumentsPattern" -> {_, _, 
   "OptionNames" -> 
    Union[First /@ Options[Gravifer`myKrylovLinearSolve], 
     First /@ Options[SparseArray`KrylovLinearSolve]]};
Gravifer`myKrylovLinearSolve[m_, b_, 
  opt : OptionsPattern[SparseArray`KrylovLinearSolve]] := 
 Module[{dim = Dimensions[m][[1]], 
    maxiter = OptionValue[MaxIterations], tol = OptionValue[Tolerance],
    x0 = OptionValue["StartingVector"], 
    resnorm = OptionValue["ResidualNormFunction"] /. Automatic -> Norm,
   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[Gravifer`myKrylovLinearSolve::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,
      If[resnorm[r] > \[Epsilon], 
       Message[Gravifer`myKrylovLinearSolve::krymit, \
\[FormalCapitalN]]]; x
   If[Precision[{m, b}] == \[Infinity], 
    Automatic | "ConjugateGradient", 
    CGfunc[dim, m, b, maxiter /. Automatic -> n, 
     tol /. Automatic -> 10^-12, 
     x0 /. Automatic -> ConstantArray[0, dim]]
   ] /; And[CheckArguments[Gravifer`myKrylovLinearSolve[m, b, opt], 2],
   (If[MatrixQ[m], True, 
     Message[Gravifer`myKrylovLinearSolve::matrix, m, 1];
     If[MatrixQ[Evaluate@m], True, 
      Message[Gravifer`myKrylovLinearSolve::krynfa, m, 1]]; False]),
   (VectorQ[b] \[Or] 
     If[MatrixQ[b], True, 
      Message[Gravifer`myKrylovLinearSolve::matrix, b, 2]; False]),
   ((OptionValue[MaxIterations] == Automatic) \[Or] 
     If[IntegerQ@OptionValue[MaxIterations], True, 
       OptionValue[MaxIterations]]; False]),
   ((OptionValue["StartingVector"] == Automatic) \[Or] 
     If[VectorQ[OptionValue["StartingVector"], NumberQ], True, 
      Message[Gravifer`myKrylovLinearSolve::kryivec]; False]),
    Automatic | "ConjugateGradient", True,
    "SteepestDescent" | "Lanczos" | "Arnoldi" | "MinRES" | "GMRES" | 
     "CGR" | "CGNR" | "CGLE" | "QMR" | "BiCG" | "BiCGStab", 
    Message[Gravifer`myKrylovLinearSolve::krynimp]; False,
    _, Message[Gravifer`myKrylovLinearSolve::kryme, 
     OptionValue["Method"]]; False

On the cases from his answer,

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} = 
 Gravifer`myKrylovLinearSolve[A, b, Method -> "ConjugateGradient", 
   "ResidualNormFunction" -> ((Sow[Norm[#]]) &)] // 
ListLinePlot[{reap[[1]], myreap[[1]]}, 
 ScalingFunctions -> "Log"]
{sol, reap} = 
  LinearSolve[N@A, N@b, 
    Method -> {"Krylov", Method -> "ConjugateGradient", 
      "ResidualNormFunction" -> ((Sow[Norm[#]]) &)}] // Reap;
{mysol, myreap} = 
 Gravifer`myKrylovLinearSolve[N@A, N@b, Method -> "ConjugateGradient",
    "ResidualNormFunction" -> ((Sow[Norm[#]]) &)] // 
ListLinePlot[{reap[[1]], myreap[[1]]}, 
 ScalingFunctions -> "Log"]

The results are

ListLinePlot[{reap[[1]], myreap[[1]]},
ScalingFunctions -> "Log"]

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

  • $\begingroup$ 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. $\endgroup$
    – user293787
    Nov 16, 2022 at 4:56
  • $\begingroup$ @user293787 I did some small scale tests, and I really assumed tha was the case. I'm a little confused right now... darn it $\endgroup$
    – Gravifer
    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;
       r = A.x-b;
       p = {r,First[p,Nothing]} // Developer`ToPackedArray;
       p[[1]] = LeastSquares[p.A.Transpose[p],-p.r].p;
       x += p[[1]];

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

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

Monitor progress:


enter image description here

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} *)

Norm[A.LinearSolve[A,b,Method->{"Krylov",Method->"ConjugateGradient"}]-b] // AbsoluteTiming
(* {11.5928, 2.43374*10^-10} *)
  • $\begingroup$ Thanks! I do have a toy CG by myself, but I indeed forgot to use AbsoluteTiming to benchmark it. $\endgroup$
    – Gravifer
    Nov 16, 2022 at 1:08

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.