Using the Krylov method for Solve: Speeding up a SparseArray calculation

Question

I'm trying to implement this Total Variation Regularized Numerical Differentiation (TVDiff) code in MMA (which I found through this SO answer): essentially I want to differentiate noisy data. The full paper behind the idea is available from LANL. For a related idea, see this Wikipedia article.

The problem I am currently having is the very long time it takes to solve the Solve function. The full TVD function is as follows:

TVD[data_, dx_] := 
 Module[{n, ep, c, DD, DDT, A, AT, u0, u, ATb, ofst, kernel, alph, ng,
    stopTol, d, iter, xs, i = 0, q, l, g, tol, maxit, p, time, s},
  n = Length@data;
  ep = 1 10^-6;
  c = Table[1/dx, {n + 1}];
  DD = SparseArray[{Band[{1, 1}] -> -c, Band[{1, 2}] -> c}, {n, 
     n + 1}];
  DDT = DD\[Transpose];

  A[list_] := dx Drop[(Accumulate[list] - 1/2 (list + list[[1]])), 1] ;
  AT[list_] := 
   dx ( Total[list] Table[1, {n + 1}] - 
      Join[{Total[list]/2}, (Accumulate[list] - list/2)]);

  ofst = data[[1]];
  ATb = AT[ofst - data];

  kernel[m_] := 
   Join[{0}, Table[Exp[-1] BesselI[n, 1], {n, -m, m}], {0}]/
    Total[Join[{0}, Table[Exp[-1] BesselI[n, 1], {n, -m, m}], {0}]];

  u0 = Join[{0}, Differences[data], {0}];
  u = SparseArray[u0];

  alph = StandardDeviation@ListConvolve[kernel[2], data]/
    StandardDeviation[data];

  ng = Infinity;
  d = 0;
  stopTol = .05;

  iter = 100;
  xs = Table[Symbol["x" <> ToString@i], {i, n + 1}];
  For[i = 0, i < iter, i++,
   q = SparseArray[Band[{1, 1}] -> 1/Sqrt[(DD.u)^2 + ep], {n, n}];
   l = dx DDT.q.DD;
   g = AT[ A[ u ] ] + ATb + alph l.u;

   tol = 10^-4;
   maxit = 1;
   (* preconditioner *)
   p = alph SparseArray[Band[{1, 1}] -> Diagonal[l], {n + 1, n + 1}];
   time = 
    AbsoluteTiming[
     s = xs /. 
        First[Solve[Thread[alph l.xs + AT[ A[ xs ] ] == g], xs, 
          Method -> {"Krylov", Method -> "ConjugateGradient", 
            "Preconditioner" -> (p.# &), Tolerance -> tol, 
            MaxIterations -> maxit}]];];
   u = u - s;
   If[Norm[s]/Norm[u] < stopTol, Break[];];
   ];
  u
  ]

Note: lower the stopTol value to ensure a better resulting derivative.

For comparison, the Matlab code (which I translated to MMA, and is available from the first link) for the "solve" portion is as follows:

s = pcg( @(v) ( alph * L * v + AT( A( v ) ) ), g, tol, maxit, P )

Here, Matlab defines the solver pcg as:

pcg(A,b,tol,maxit,M) and pcg(A,b,tol,maxit,M1,M2) use symmetric positive definite preconditioner M or M = M1*M2 and effectively solve the system inv(M)*A*x = inv(M)*b for x. If M is [] then pcg applies no preconditioner. M can be a function handle mfun such that mfun(x) returns M\x.

Note also that the @(v) is a 'pure function' in Matlab terms, and is allowed as per:

A can be a function handle afun such that afun(x) returns A*x.

When I run the two codes, Matlab ends up being ~5-20 times faster than the corresponding MMA code. My MMA implementation of it uses more or less the entire cpu time on the Solve function.

I tried to find the best corresponding MMA Solver routine that matched the Matlab description via the docs and two different MathGroup archived messages. None of the options (whether given with or without quotes) seem to help at all.

For testing purposes, here is some data:

data = {4699.1`, 4728.3`, 4753.3`, 4787.4`, 4794.8`, 4817.5`, 4842.7`,
   4877.2`, 4888.2`, 4916.1`, 4933.7`, 4951.5`, 4984.1`, 4984.2`, 
  5004.`, 5031.`, 5048.1`, 5062.3`, 5083.2`, 5096.`, 5108.5`, 5140.`, 
  5142.8`, 5142.7`, 5169.1`, 5168.6`, 5165.`, 5191.8`, 5193.7`, 
  5199.4`, 5189.3`, 5213.6`, 5209.1`, 5208.5`, 5197.`, 5201.2`, 
  5184.2`, 5191.2`, 5183.7`, 5181.3`, 5183.2`, 5175.6`, 5089.9`, 
  5068.1`, 5053.9`, 5056.7`, 5063.6`, 5038.2`, 5023.9`, 5027.4`, 
  4998.8`, 4980.9`, 4961.9`, 4939.3`, 4933.`, 4897.7`, 4879.`, 4874.`,
   4857.3`, 4819.2`, 4801.6`, 4775.5`, 4754.9`, 4712.2`, 4708.3`, 
  4675.8`, 4637.1`, 4634.1`, 4582.6`, 4558.3`, 4531.`, 4507.9`, 
  4470.4`, 4445.7`, 4435.`, 4404.3`, 4383.5`, 4363.7`}

You can make it bigger (which also coincidentally crashes my MMA if too big...) by:

data = Join[data, Reverse@data, data, Reverse@data, data, 
   Reverse@data, data];

(rinse & repeat as necessary). The data looks like:

Mathematica graphics

And the function TVD[data,1/Length@data] looks like:

Mathematica graphics

So, how can I speed the solver up? Is Matlab just that much better at 'sparse array' type calculations? Did I not define the right SparseArrays? Is there a way to use LinearSolve when the matrix equation is not a simple A.x on the left hand side?

Any and all speed improvements would be great!

I haven't looked into this method and don't have time right now so I can't answer your question at the moment. However if you just want to differentiate noisy data and don't care about the approach, you can perhaps use the Savitzky-Golay method. I posted some code on MathGroup here.
You asked: <"Is Matlab just that much better at 'sparse array' type calculations? "> I do not think so. I believe that Both Mathematica's LinearSolve and Matlab's "\" end up using the UMFPACK library. cise.ufl.edu/research/sparse/umfpack and may be also SPARSKIT. I find Matematica's LinearSolve very fast. I use both Matlab and Mathematica all the time.They both use the same core linear algebra libraries below all these API's (such as intel's MKL). You just need to setup your problem so that you use the right solver and use the right form to get most of the benefits and the speed.
@NasserM.Abbasi I asked that mostly rhetorically. While matrix things in matlab are generally easier to do than in MMA, I much prefer MMA's syntax / computational abilities. I'm really hoping that I'm just missing the right options and / or understanding of how to make my problem work with MMA's Solve / LinearSolve
@tkott, when I run your code I get messages :-( SparseArray::bndtr:

Oleksandr R. · Accepted Answer · 2012-04-21 00:10:48Z

I found a way to dramatically improve the performance of this algorithm by using the undocumented function SparseArray`KrylovLinearSolve. The key advantage of this function is that it seems to be a near-analog of MATLAB's pcg, and as such accepts as a first argument either:

a square matrix, or a function generating a vector of length equal to the length of the second argument.

One may discover this by giving incorrect arguments and noting the advice given in the messages produced as a result, in much the same way as one discovers the correct arguments for any undocumented function. In this case the message is SparseArray`KrylovLinearSolve::krynfa.

You only need to change one line in your code to use it, namely:

s = SparseArray`KrylovLinearSolve[
     alph l.# + AT[A[#]] &, g, 
     Method -> "ConjugateGradient", "Preconditioner" -> (p.# &), 
     Tolerance -> tol, MaxIterations -> maxit
    ];

where maxit should preferably be Automatic (meaning 10 times the size of the system to be solved) or larger. With the data given in your question it takes a few hundred iterations to converge to a tolerance of $10^{-4}$, but each iteration is quite fast, so it seems to make more sense to adjust the tolerance than the number of iterations if performance is still an issue. However, while I didn't investigate this, needing this many iterations to converge to a relatively loose tolerance may of course be symptomatic of a poorly conditioned system, so using a different preconditioner or the biconjugate gradient stabilized method ("BiCGSTAB") could perhaps reduce the number of iterations required.

You will note that the options are exactly the same as for LinearSolve's "Krylov" method, so we may surmise that this function is probably called more or less directly by LinearSolve when Method -> "Krylov" is specified. In fact, if we assume that this is indeed the case and try

s = LinearSolve[
     alph l.# + AT[A[#]] &, g, 
     Method -> {"Krylov",
       Method -> "ConjugateGradient", "Preconditioner" -> (p.# &), 
       Tolerance -> tol, MaxIterations -> maxit
      }
    ];

we find that it works equally well, so evidently LinearSolve does in fact provide just the same functionality as pcg as far as the first argument is concerned, but without this actually being documented anywhere as far as I can tell. So, the overall conclusion is that you can just use LinearSolve directly after all.

Stupid undocumented functions. Thanks Oleksandr! I'll play around with this on Monday :)
Finally got around to checking this darned thing. Your solution of using LinearSolve + the pure function worked like a charm! Thanks

Nasser · Answer 2 · 2012-04-18 17:39:10Z

This is not an answer, but may be it is :) as I do not have time to fully understand the question, just picked up few terms, but just in case, I thought I mention this.

Mathematica 8 already has Krylov method in LinearSolve ! so, if you are just looking to use these methods to solve Ax=b, it is already there.

Here is an example from my code (I used these in a demo for solving poisson 2D pde)

Which[preconditioner == "ILU0",
  x = LinearSolve[A, rightHandVector, 
    Method -> {"Krylov", Method -> nonStationarySolver, 
      "Preconditioner" -> preconditioner, MaxIterations -> Automatic, 
      Tolerance -> 10^-toleranceConstant}],
  True,
  x = LinearSolve[A, rightHandVector, 
    Method -> {"Krylov", Method -> nonStationarySolver, 
      "Preconditioner" -> {preconditioner, "FillIn" -> fillIn}, 
      MaxIterations -> Automatic, Tolerance -> 10^-toleranceConstant}]
  ];

In the above A is sparse 2D matrix and rightHandVector is the 'b' vector in Ax=b

There are many other submethods for Krylov and options If this is what you are asking for, then you can check my demo for more example, I also implement Preconditioned Conjugate Gradients Method (what Matlab pcg does).

If this is NOT what you are asking for, I can delete this. (the comment was too small to write all this in).

ps. the code for this is here http://12000.org/my_notes/mma_demos/poisson_2D/index.htm

I would love to use LinearSolve. however, my left hand side "matrix" is really: l.x + AT[ A[ x ] ] (i.e., it is not of a form M.x==b, but in some sense f[M]==b, where the functions f takes care of M.x. If there is a way to do make this work with MMA's LinearSolve, I'm all for it.
Simple question: where did you find the list of methods available for LinearSolve?
@rcollyer, Oliver Ruebenkoenig from WRI was kind enough to reply to a post I had at Mathgroup. If you google "A function to do incomplete LU decomposition with a drop tolerance" and select the second or third link that comes up, you'll see the thread and his reply. I know some of these are hard to find out from the current documentation center.
Difficult?!? I'd say it borders on impossible. Maybe it is in need of a question?
@rcollyer Actually I link to that discussion in my question :) (there are two links right before "MathGroup")

asked	1 year ago
viewed	460 times
active	1 year ago

Using the Krylov method for Solve: Speeding up a SparseArray calculation

2 Answers

Your Answer

Not the answer you're looking for? Browse other questions tagged performance-tuning equation-solving matrix sparse-arrays or ask your own question.

Community Bulletin

Linked

Using the Krylov method for Solve: Speeding up a SparseArray calculation

2 Answers

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged performance-tuning equation-solving matrix sparse-arrays or ask your own question.

Community Bulletin

Linked

Related