2

From the documentation, we learn that:

LeastSquares[m,b], when b is a vector, is equivalent to ArgMin[Norm[m.x-b],x]. I am wondering whether, at the level of performances (say, at the level of convergence and speed), is there a reason to prefer one method or the other when we have to deal with large and sparse matrices and vectors (and we are solving the problem numerically). I know that the question is a bit general, but any clue is more than welcome.

UPDATE

Having decided to use LeastSquare I have a further question: is LeastSquare parallelizable? For example, by doing a compiled version of the function, with

CompilationTarget -> "C" , RuntimeAttributes -> {Listable} , Parallelization -> True

Should I expect some significative speed up?

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
Dario Rosa
  • 355
  • 1
  • 9
  • 4
    LeastSquares is very much preferred for large numerical matrices (both sparse and full). It's specifically optimized that. – Sjoerd Smit Jun 20 '20 at 10:08
  • Thanks! just noather point: Do you know how I can force LeastSquares to look for real solutions only? – Dario Rosa Jun 20 '20 at 11:59
  • 4
    I suspect LeastSquares is just solving the normal equations (I think that's the term), which is the more efficient version of multiplying by the pseudoinverse matrix. So there is no way to place a restriction on it. – Daniel Lichtblau Jun 20 '20 at 14:12
  • I see. @DanielLichtblau, thank you – Dario Rosa Jun 20 '20 at 14:17
  • 2
    @DanielLichtblau Oh, I very much hope that LeastSquares avoids the normal equations; those typically have a much larger condition number. More stable alorithms are based on QR-decomposition or SVD. – Henrik Schumacher Jun 20 '20 at 14:17
  • 3
    @HenrikSchumacher It's not solving them directly. I believe it is actually using SVD, which is slower than QR but also safer. – Daniel Lichtblau Jun 20 '20 at 20:46

1 Answers1

4

LeastSquare seems to use MKL lib internally and is already parallelized. Run the following code and observe cpu load:

SeedRandom[1] ;
size = 3000 ;
m = RandomReal[{-1,1},{size,size}] ;
b = RandomReal[{-1,1},size] ;
LeastSquares[m,b] ; // AbsoluteTiming

Check your parallel options and compare with observed cpu load:

SystemOptions["ParallelOptions"]

You will get no speed up from Compile if you compiling only LeastSquares function, it will not be compiled, there will be a call to MainEvaluate:

<<CompiledFunctionTools`
CompilePrint[Compile[{{m,_Real,2},{b,_Real,1}},LeastSquares[m,b]]]

To get a potential speed up, all functions inside Compile should be compilable. If you have a function like Compile[{...},code1,LeastSquares[...],code2], you can try to compile code1 and code2 separately.

You can try low level functions to replace LeastSquares, see LowLevelLinearAlgebra/guide/BLASGuide in docs. This might be helpful if your matrix have special structure, e.g. triangular. More low level function:

?LinearAlgebra`*`*
I.M.
  • 2,926
  • 1
  • 13
  • 18
  • Nice answer! I will try the low level functions you suggest! – Dario Rosa Jun 21 '20 at 04:49
  • 1
    @DarioRosa, thanks, also note, both options RuntimeAttributes -> {Listable}, Parallelization -> True should present if your want speed up from Compile parallelization and compiled function invocation is in it's listable form. This is similar to do-loop parallelization in openMP. – I.M. Jun 21 '20 at 04:58