5

I'm converting a piece of matlab code to python but I am having difficulty recreating the solution it finds for an overdetermined sparse matrix. The original code used the overloaded left division operator but I have recreated the results using linsolve(). I have compared the result with several scikit solvers (OLS, ridge regression, lasso, etc) but have not managed to match it. Lasso came closest as the matlab coefficients are sparse so I think it must be something like that.

I have debugged it and used type to try and access the source code, with no success. Does anyone know what solver it uses or how I could get a peak "under the hood"?

jonathanbyrn
  • 151
  • 3

2 Answers2

5

MATLAB's linsolve() function uses QR factorization with column pivoting to find a least squares solution to your over determined problem.

You can use

[x,r]=linsolve(A,b)

to get an estimate of the rank of $A$. I believe that if $A$ is rank deficient, then MATLAB will set all of the free variables to $0$ in computing the least squares solution. If this is happening, it could explain why MATLAB's least squares solution is sparse.

If your matrix is of full column rank, then the least squares solution should be unique and your Python code should be finding essentially the same solution.

If your matrix not of full column rank, then there will be infinitely many least squares solutions, and it would be no surprise that your python code found a very different solution. As long as $\| Ax-b \|_{2}$ is the same for the solution obtained by your code and by MATLAB, you've got no real reason to complain.

The most complicated situation occurs when your matrix is badly conditioned so that it is on the very edge of not having full column rank. If you're in this situation then there may be many solutions that have very nearly equal values of $\| Ax -b \|_{2}$ and very different values of $x$.

I would suggest three things to look at:

  1. Compare the values of $\| Ax-b \|_{2}$ for MATLAB's solution and the solutions produced by your Python code. They should be essentially equal.

  2. Use [x,r]=linsolve(A,b) in MATLAB to get a least squares solution together with an estimate of the rank of $A$. This will give you some indication of whether $A$ is rank deficient.

  3. Take a reasonably small test problem and compute the condition number of your $A$ matrix (in MATLAB, cond(A) will do this.) This should tell you whether you have a matrix that is effectively of full column rank or not.

The answers to these questions might help us in directing you further.

Brian Borchers
  • 18,719
  • 1
  • 39
  • 70
  • Also see the documentation for linsolve(), which discusses what the function is doing. – Brian Borchers Jul 13 '13 at 01:43
  • For a fairly complete discussion of the Q-less QR factorization approach that MATLAB uses for sparse overdetermined least squares problems, see http://www.mathworks.com/help/matlab/math/sparse-matrix-operations.html – Brian Borchers Jul 13 '13 at 05:12
  • By the way, at least in MATLAB R2012a (we're not running a newer version), linsolve() only takes dense matrices as inputs, although backslash uses the Q-less QR approach. – Brian Borchers Jul 13 '13 at 05:28
  • Thank you so much for your help. Just to clarify I am adapting echo state network code. I am using ∥Ax−b∥2 to compare the different methods and I would be happy enough if they were close but unfortunately they are not. A is indeed rank deficient (rank = 11) so that explains the sparse result. when I find the condition number of A it is very large( 2.5290e+17) so this means it is not a well conditioned matrix. – jonathanbyrn Jul 17 '13 at 14:11
  • That clarifies the problem a bit. I'm still not clear on whether your MATLAB code is using linsolve() on $A$ stored in dense format, or whether you're using backslash on $A$ stored in sparse format. Slightly different algorithms are used in these two cases. How does "whos" in MATLAB describe your $A$ matrix? – Brian Borchers Jul 17 '13 at 14:35
3

The answer https://scicomp.stackexchange.com/a/1459/276 provides some useful advice.

If you insert the line spparms('spumoni',1); before your sparse solve, it will output the decisions made by the algorithm. Chances are, it uses a sparse QR factorization.

Community
  • 1
Geoff Oxberry
  • 30,394
  • 9
  • 64
  • 127
  • The example from that thread worked for me but when I tried it on my own matrix I saw nothing. I think it might be because the matrix is rank deficient. I will try out a QR factorization technique and see if I get the same results, thanks for your help! – jonathanbyrn Jul 16 '13 at 16:33
  • Did you try reverting the code to use backslash instead of linsolve, as you mentioned you had tried earlier in your question? – Geoff Oxberry Jul 16 '13 at 17:32