3

I am trying to solve a system of linear equations in the least-squares-style. Using armadillo and its solve function I want to calculate the three coefficients of a parabolic fit.

vec coeffs = solve(CtC, Ctb)

with CtC=

1.0e+009 *       
+--------------------------------+
|     2.0878    0.0221    0.0002 |
|     0.0221    0.0002    0.0000 |
|     0.0002    0.0000    0.0000 |
+--------------------------------+

and Ctb=

+------------+
|    -0.6163 |
|    -0.0065 |
|    -0.0001 |
+------------+

Apparently solve() does not manage to solve it, even Matlab warns:

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate.
RCOND = 1.303968e-022.

Is there any workaround or a more powerful/sophisticated method in armadillo or c++ in general? Thanks

TemplateRex
  • 69,038
  • 19
  • 164
  • 304
bogus
  • 867
  • 6
  • 14
  • 24

2 Answers2

2

You cannot invert a matrix with determinant equal to zero. In your case, the determinant is equal to -8e-12, which is not actually equal to zero but which gives as inverse matrix (using R)

1.0e-009 *
              [,1]          [,2]    [,3]
[1,] -3.815763e-16  6.806617e-14    5000
[2,]  5.111864e-14  5.000000e+03 -552500
[3,]  5.000000e+03 -5.525000e+05 8856250

Note that this inverse matrix has both very large and very small numbers. This explains your warning about the "close to singular" and "badly scaled". Depending on the underlying floating point library, you might get very bad rounding errors.

You could go back to your application and see if there are any variables representing a column or row that are many orders of magnitude apart from each other. If so, you could rescale those. E.g. if there is a factor of 1 million between them, and they represent length, express one in kilometers and the other in millimeters.

TemplateRex
  • 69,038
  • 19
  • 164
  • 304
  • Yes, Im aware of that, it was to be expected according to [this](http://arma.sourceforge.net/docs.html). I was thinking of an implementation to compute the coefficients using a pseudoinverse such as Moore-Penrose. Does that an implementation like that exist? – bogus Jan 23 '13 at 13:44
  • @bogus from your notation, I gather that your original problem is `C x = b` and that `C` is not a square matrix. This means that you were trying to solve `CtC x = Ct b`, i.e. you multiply both sides of your original equation by `Ct` (C tranpose). This is already a generalized inverse, so you won't gain anything. – TemplateRex Jan 23 '13 at 13:48
  • Exactly, C is not squadratic, that why I called it "least squares style". I found a well-working solution for my practical application, which makes sense in my context. Thanks anyways! – bogus Jan 23 '13 at 13:56
  • @bogus So what did you do? This might help others who are having the same problem... – Martin B Jan 23 '13 at 14:09
  • The determinant and the inverse are not the right ones, you forgot the 1e9 factor in front of the matrix. However, this does not change the problem and your answer is still valid. – Dr_Sam Jan 23 '13 at 15:22
2

The problem that you are facing is that the condition number of your matrix is very high, which might make your results meaningless due to rounding errors.

As proposed by rhalbersma, you can try to better scale the different entries to reduce this problem. If you have no such solution, then you have to use a preconditioner. There are a lot of them (some very complicated) so you should find one that fits your problem.

However, for this particular system that you have at hand, I would simply invert the equations (1) and (3), which would yield a lower triangular system, that is easy to solve for Matlab (no need for iterative methods which can be very sensitive to the condition number. In your case, since the matrix is symmetric, Matlab might try to solve it using cg, which is an iterative method).

By the way, is your problem known to yield near-singular matrices? If not, I would check for a problem (or bug) before the assembly of the matrix (the matrix might be non-singular thanks to rounding errors).

Edit: On my matlab 2009, I do not get any warning while trying to solve this system (using the \ command). However, if I try to use the command pcg, I get them!

Dr_Sam
  • 1,818
  • 18
  • 24
  • +1 for mentioning preconditioners. I had completely forgotten about those (it's 10+ years since my numerical analysis courses) – TemplateRex Jan 23 '13 at 15:25