numpy (matrix solver) - python vs. matlab

[someone]

Hi,

Notice cross-post, I hope you bear over with me for doing that (and I
imagine that some of you also like python in the matlab-group like
myself)...

------------------------------------------
Python vs. Matlab:
------------------------------------------

Python:
========
from numpy import matrix
from numpy import linalg
A = matrix( [[1,2,3],[11,12,13],[21,22,23]] )
print "A="
print A
print "A.I (inverse of A)="
print A.I

A.I (inverse of A)=
[[ 2.81466387e+14 -5.62932774e+14 2.81466387e+14]
[ -5.62932774e+14 1.12586555e+15 -5.62932774e+14]
[ 2.81466387e+14 -5.62932774e+14 2.81466387e+14]]


Matlab:
========

A=[1 2 3; 11 12 13; 21 22 23]

A =

1 2 3
11 12 13
21 22 23

inv(A)

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

ans =

1.0e+15 *

0.3002 -0.6005 0.3002
-0.6005 1.2010 -0.6005
0.3002 -0.6005 0.3002

------------------------------------------
Python vs. Matlab:
------------------------------------------

So Matlab at least warns about "Matrix is close to singular or badly
scaled", which python (and I guess most other languages) does not...

Which is the most accurate/best, even for such a bad matrix? Is it
possible to say something about that? Looks like python has a lot more
digits but maybe that's just a random result... I mean.... Element 1,1 =
2.81e14 in Python, but something like 3e14 in Matlab and so forth -
there's a small difference in the results...

With python, I would also kindly ask about how to avoid this problem in
the future, I mean, this maybe means that I have to check the condition
number at all times before doing anything at all ? How to do that?

I hope you matlabticians like this topic, at least I myself find it
interesting and many of you probably also program in some other language
and then maybe you'll find this worthwhile to read about.

[Kiuhnm]

Post by someone
Hi,
Notice cross-post, I hope you bear over with me for doing that (and I
imagine that some of you also like python in the matlab-group like
myself)...
------------------------------------------
------------------------------------------
========
from numpy import matrix
from numpy import linalg
A = matrix( [[1,2,3],[11,12,13],[21,22,23]] )
print "A="
print A
print "A.I (inverse of A)="
print A.I
A.I (inverse of A)=
[[ 2.81466387e+14 -5.62932774e+14 2.81466387e+14]
[ -5.62932774e+14 1.12586555e+15 -5.62932774e+14]
[ 2.81466387e+14 -5.62932774e+14 2.81466387e+14]]
========

A=[1 2 3; 11 12 13; 21 22 23]

A =
1 2 3
11 12 13
21 22 23

inv(A)

Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 1.067522e-17.
ans =
1.0e+15 *
0.3002 -0.6005 0.3002
-0.6005 1.2010 -0.6005
0.3002 -0.6005 0.3002
------------------------------------------
------------------------------------------
So Matlab at least warns about "Matrix is close to singular or badly
scaled", which python (and I guess most other languages) does not...

A is not just close to singular: it's singular!

Post by someone
Which is the most accurate/best, even for such a bad matrix? Is it
possible to say something about that? Looks like python has a lot more
digits but maybe that's just a random result... I mean.... Element 1,1 =
2.81e14 in Python, but something like 3e14 in Matlab and so forth -
there's a small difference in the results...

Both results are *wrong*: no inverse exists.

Post by someone
With python, I would also kindly ask about how to avoid this problem in
the future, I mean, this maybe means that I have to check the condition
number at all times before doing anything at all ? How to do that?

If cond(A) is high, you're trying to solve your problem the wrong way.
You should try to avoid matrix inversion altogether if that's the case.
For instance you shouldn't invert a matrix just to solve a linear system.

Kiuhnm

[someone]

Post by Kiuhnm

Post by someone
So Matlab at least warns about "Matrix is close to singular or badly
scaled", which python (and I guess most other languages) does not...

A is not just close to singular: it's singular!

Ok. When do you define it to be singular, btw?

Post by Kiuhnm

Post by someone
Which is the most accurate/best, even for such a bad matrix? Is it
possible to say something about that? Looks like python has a lot more
digits but maybe that's just a random result... I mean.... Element 1,1 =
2.81e14 in Python, but something like 3e14 in Matlab and so forth -
there's a small difference in the results...

Both results are *wrong*: no inverse exists.

What's the best solution of the two wrong ones? Best least-squares
solution or whatever?

Post by Kiuhnm

Post by someone
With python, I would also kindly ask about how to avoid this problem in
the future, I mean, this maybe means that I have to check the condition
number at all times before doing anything at all ? How to do that?

If cond(A) is high, you're trying to solve your problem the wrong way.

So you're saying that in another language (python) I should check the
condition number, before solving anything?

Post by Kiuhnm
You should try to avoid matrix inversion altogether if that's the case.
For instance you shouldn't invert a matrix just to solve a linear system.

What then?

Cramer's rule?

[Nasser M. Abbasi]

Post by someone
Ok. When do you define it to be singular, btw?

There are things you can see right away about a matrix A being singular
without doing any computation. By just looking at it.

For example, If you see a column (or row) being a linear combination of
other column(s) (or row(s)) then this is a no no.

In your case you have

1 2 3
11 12 13
21 22 23

You can see right away that if you multiply the second row by 2, and
subtract from that one times the first row, then you obtain the third row.

Hence the third row is a linear combination of the first row and the
second row. no good.

When you get a row (or a column) being a linear combination of others
rows (or columns), then this means the matrix is singular.

--Nasser

[Kiuhnm]

Post by someone

Post by Kiuhnm

Post by someone
So Matlab at least warns about "Matrix is close to singular or badly
scaled", which python (and I guess most other languages) does not...

A is not just close to singular: it's singular!

Ok. When do you define it to be singular, btw?

Post by Kiuhnm

Post by someone
Which is the most accurate/best, even for such a bad matrix? Is it
possible to say something about that? Looks like python has a lot more
digits but maybe that's just a random result... I mean.... Element 1,1 =
2.81e14 in Python, but something like 3e14 in Matlab and so forth -
there's a small difference in the results...

Both results are *wrong*: no inverse exists.

What's the best solution of the two wrong ones? Best least-squares
solution or whatever?

Trust me. They're both so wrong that it doesn't matter.
Have a look at A*inv(A) and inv(A)*A and you'll see by yourself.

Post by someone

Post by Kiuhnm

Post by someone
With python, I would also kindly ask about how to avoid this problem in
the future, I mean, this maybe means that I have to check the condition
number at all times before doing anything at all ? How to do that?

If cond(A) is high, you're trying to solve your problem the wrong way.

So you're saying that in another language (python) I should check the
condition number, before solving anything?

Yes, unless you already know that it will always be low by design.

Post by someone

Post by Kiuhnm
You should try to avoid matrix inversion altogether if that's the case.
For instance you shouldn't invert a matrix just to solve a linear system.

What then?

Look at the documentation of the library you're using.

Post by someone
Cramer's rule?

Surprisingly, yes. That's an option. See
"A condensation-based application of Cramerʼs rule for solving
large-scale linear systems"
Popular linear codes are based on Gaussian elimination or some iterative
method, though.

Kiuhnm

[Russ P.]

Post by someone

Post by Kiuhnm

Post by someone
So Matlab at least warns about "Matrix is close to singular or badly
scaled", which python (and I guess most other languages) does not...

A is not just close to singular: it's singular!

Ok. When do you define it to be singular, btw?

Post by Kiuhnm

Post by someone
Which is the most accurate/best, even for such a bad matrix? Is it
possible to say something about that? Looks like python has a lot more
digits but maybe that's just a random result... I mean.... Element 1,1 =
2.81e14 in Python, but something like 3e14 in Matlab and so forth -
there's a small difference in the results...

Both results are *wrong*: no inverse exists.

What's the best solution of the two wrong ones? Best least-squares
solution or whatever?

Post by Kiuhnm

Post by someone
With python, I would also kindly ask about how to avoid this problem in
the future, I mean, this maybe means that I have to check the condition
number at all times before doing anything at all ? How to do that?

If cond(A) is high, you're trying to solve your problem the wrong way.

So you're saying that in another language (python) I should check the
condition number, before solving anything?

Post by Kiuhnm
You should try to avoid matrix inversion altogether if that's the case.
For instance you shouldn't invert a matrix just to solve a linear system.

What then?
Cramer's rule?

If you really want to know just about everything there is to know
about a matrix, take a look at its Singular Value Decomposition (SVD).
I've never used numpy, but I assume it can compute an SVD.

[Eelco]

There is linalg.pinv, which computes a pseudoinverse based on SVD that
works on all matrices, regardless of the rank of the matrix. It merely
approximates A*A.I = I as well as A permits though, rather than being
a true inverse, which may not exist.

Anyway, there are no general answers for this kind of thing. In all
non-textbook problems I can think of, the properties of your matrix
are highly constrained by the problem you are working on; which
additional tests are required to check for corner cases thus depends
on the problem. Often, if you have found an elegant solution to your
problem, no such corner cases exist. In that case, MATLAB is just
wasting your time with its automated checks.

[someone]

Post by Russ P.

Post by someone

Post by Kiuhnm
You should try to avoid matrix inversion altogether if that's the case.
For instance you shouldn't invert a matrix just to solve a linear system.

What then?
Cramer's rule?

If you really want to know just about everything there is to know
about a matrix, take a look at its Singular Value Decomposition (SVD).

I know a bit about SVD - I used it for a short period of time in Matlab,
though I'm definately not an expert in it and I don't understand the
whole theory with orthogality behind making it work so elegant as it
is/does work out.

Post by Russ P.
I've never used numpy, but I assume it can compute an SVD.

I'm making my first steps now with numpy, so there's a lot I don't know
and haven't tried with numpy...

[Nasser M. Abbasi]

Post by someone
I would also kindly ask about how to avoid this problem in
the future, I mean, this maybe means that I have to check the condition
number at all times before doing anything at all ? How to do that?

I hope you'll check the condition number all the time.

You could be designing a building where people will live in it.

If do not check the condition number, you'll end up with a building that
will fall down when a small wind hits it and many people will die all
because you did not bother to check the condition number when you solved
the equations you used in your design.

Also, as was said, do not use INV(A) directly to solve equations.

--Nasser

[someone]

Post by Nasser M. Abbasi

Post by someone
I would also kindly ask about how to avoid this problem in
the future, I mean, this maybe means that I have to check the condition
number at all times before doing anything at all ? How to do that?

I hope you'll check the condition number all the time.

So how big can it (cond-number) be before I should do something else?
And what to do then? Cramers rule or pseudoinverse or something else?

Post by Nasser M. Abbasi
You could be designing a building where people will live in it.
If do not check the condition number, you'll end up with a building that
will fall down when a small wind hits it and many people will die all
because you did not bother to check the condition number when you solved
the equations you used in your design.
Also, as was said, do not use INV(A) directly to solve equations.

In Matlab I used x=A\b.

I used inv(A) in python. Should I use some kind of pseudo-inverse or
what do you suggest?

[Nasser M. Abbasi]

Post by someone

Post by Nasser M. Abbasi
Also, as was said, do not use INV(A) directly to solve equations.

In Matlab I used x=A\b.

good.

Post by someone
I used inv(A) in python. Should I use some kind of pseudo-inverse or
what do you suggest?

I do not use python much myself, but a quick google showed that pyhton
scipy has API for linalg, so use, which is from the documentation, the
following code example

X = scipy.linalg.solve(A, B)

But you still need to check the cond(). If it is too large, not good.
How large and all that, depends on the problem itself. But the rule of
thumb, the lower the better. Less than 100 can be good in general, but I
really can't give you a fixed number to use, as I am not an expert in
this subjects, others who know more about it might have better
recommendations.

--Nasser

[Kiuhnm]

Post by Nasser M. Abbasi
But you still need to check the cond(). If it is too large, not good.
How large and all that, depends on the problem itself. But the rule of
thumb, the lower the better. Less than 100 can be good in general, but I
really can't give you a fixed number to use, as I am not an expert in
this subjects, others who know more about it might have better
recommendations.

Alas, there's no fixed number and as if that wasn't enough, there are
many condition numbers, each one with different properties. For
instance, the Skeel condition number is scale-invariant and it's useful
when a matrix is ill-conditioned just because its rows are out of scale.

Kiuhnm

[someone]

Post by Nasser M. Abbasi
I do not use python much myself, but a quick google showed that pyhton
scipy has API for linalg, so use, which is from the documentation, the
following code example
X = scipy.linalg.solve(A, B)
But you still need to check the cond(). If it is too large, not good.
How large and all that, depends on the problem itself. But the rule of
thumb, the lower the better. Less than 100 can be good in general, but I
really can't give you a fixed number to use, as I am not an expert in
this subjects, others who know more about it might have better
recommendations.

Ok, that's a number...

Anyone wants to participate and do I hear something better than "less
than 100 can be good in general" ?

If I don't hear anything better, the limit is now 100...

What's the limit in matlab (on the condition number of the matrices), by
the way, before it comes up with a warning ???

[Russ P.]

Post by someone

Post by Nasser M. Abbasi
I do not use python much myself, but a quick google showed that pyhton
scipy has API for linalg, so use, which is from the documentation, the
following code example
X = scipy.linalg.solve(A, B)
But you still need to check the cond(). If it is too large, not good.
How large and all that, depends on the problem itself. But the rule of
thumb, the lower the better. Less than 100 can be good in general, but I
really can't give you a fixed number to use, as I am not an expert in
this subjects, others who know more about it might have better
recommendations.

Ok, that's a number...
Anyone wants to participate and do I hear something better than "less
than 100 can be good in general" ?
If I don't hear anything better, the limit is now 100...
What's the limit in matlab (on the condition number of the matrices), by
the way, before it comes up with a warning ???

The threshold of acceptability really depends on the problem you are
trying to solve. I haven't solved linear equations for a long time,
but off hand, I would say that a condition number over 10 is
questionable.

A high condition number suggests that the selection of independent
variables for the linear function you are trying to fit is not quite
right. For a poorly conditioned matrix, your modeling function will be
very sensitive to measurement noise and other sources of error, if
applicable. If the condition number is 100, then any input on one
particular axis gets magnified 100 times more than other inputs.
Unless your inputs are very precise, that is probably not what you
want.

Or something like that.

[someone]

Post by Russ P.

Post by someone
What's the limit in matlab (on the condition number of the matrices), by
the way, before it comes up with a warning ???

The threshold of acceptability really depends on the problem you are
trying to solve. I haven't solved linear equations for a long time,
but off hand, I would say that a condition number over 10 is
questionable.

Anyone knows the threshold for Matlab for warning when solving x=A\b ? I
tried "edit slash" but this seems to be internal so I cannot see what
criteria the warning is based upon...

Post by Russ P.
A high condition number suggests that the selection of independent
variables for the linear function you are trying to fit is not quite
right. For a poorly conditioned matrix, your modeling function will be
very sensitive to measurement noise and other sources of error, if
applicable. If the condition number is 100, then any input on one
particular axis gets magnified 100 times more than other inputs.
Unless your inputs are very precise, that is probably not what you
want.
Or something like that.

Ok. So it's like a frequency-response-function, output divided by input...

[Eric Diaz]

Regarding the condition number (K). I'm not an expert but I did put together a presentation with a slide on it. One way of looking at it is how sensitive parameter solutions are to noise. Indeed, when K is high then even the tiniest bit of noise can cause huge changes in parameter estimates. This is, in essence, the gain of sensitivity to relative error in the signal of our fitting model. From my past investigations, I believe that one way to get a practical sense for the impact of your condition number on your parameter estimates is to use the log base 2 of K. This result tells you the "digits of lost precision". So, with a K=1, there is no loss of precision, with a K~50 there is loss of ~5 digits of precision, and with a K~1000 there is loss of ~10 digits of precision. Obviously without precision, there is no reliability in the repeatability of your solution.

[Steven_Lord]

Post by Russ P.

Post by someone

Post by Nasser M. Abbasi
I do not use python much myself, but a quick google showed that pyhton
scipy has API for linalg, so use, which is from the documentation, the
following code example
X = scipy.linalg.solve(A, B)
But you still need to check the cond(). If it is too large, not good.
How large and all that, depends on the problem itself. But the rule of
thumb, the lower the better. Less than 100 can be good in general, but I
really can't give you a fixed number to use, as I am not an expert in
this subjects, others who know more about it might have better
recommendations.

Ok, that's a number...
Anyone wants to participate and do I hear something better than "less
than 100 can be good in general" ?
If I don't hear anything better, the limit is now 100...
What's the limit in matlab (on the condition number of the matrices), by
the way, before it comes up with a warning ???

I'm not going to answer, and the reason why is that saying "the limit is X"
may lead you to believe that "as long as my condition number is less than X,
I'm safe." That's not the case. The threshold above which MATLAB warns is
the fence that separates the tourists from the edge of the cliff of
singularity -- just because you haven't crossed the fence doesn't mean
you're not in danger of tripping and falling into the ravine.

Post by Russ P.
The threshold of acceptability really depends on the problem you are
trying to solve.

I agree with this statement, although you generally don't want it to be too
big. The definition of "too big" is somewhat fluid, though.

Post by Russ P.
I haven't solved linear equations for a long time,
but off hand, I would say that a condition number over 10 is
questionable.

That seems pretty low as a general bound IMO.

Post by Russ P.
A high condition number suggests that the selection of independent
variables for the linear function you are trying to fit is not quite
right. For a poorly conditioned matrix, your modeling function will be
very sensitive to measurement noise and other sources of error, if
applicable. If the condition number is 100, then any input on one
particular axis gets magnified 100 times more than other inputs.
Unless your inputs are very precise, that is probably not what you
want.
Or something like that.

Russ, you and the OP (and others) may be interested in one of the books that
Cleve Moler has written and made freely available on the MathWorks website:

http://www.mathworks.com/moler/

The chapter Linear Equations in "Numerical Computing with MATLAB" includes a
section (section 2.9, starting on page 14 if I remember correctly) that
discusses norm and condition number and gives a more formal statement of
what you described. The code included in the book is written in MATLAB, but
even if you don't use MATLAB (since I know this has been cross-posted to
comp.lang.python) there's plenty of good, crunchy mathematics in that
section.

[someone]

Post by Steven_Lord
Russ, you and the OP (and others) may be interested in one of the books
that Cleve Moler has written and made freely available on the MathWorks
http://www.mathworks.com/moler/
The chapter Linear Equations in "Numerical Computing with MATLAB"
includes a section (section 2.9, starting on page 14 if I remember
correctly) that discusses norm and condition number and gives a more
formal statement of what you described. The code included in the book is
written in MATLAB, but even if you don't use MATLAB (since I know this
has been cross-posted to comp.lang.python) there's plenty of good,
crunchy mathematics in that section.

I use matlab more than python. I just want to learn python, which seems
extremely powerful and easy but yet, I'm a python beginner.

Thank you very much for the reference, Mr. Lord. I'll look closely at
that Moler-book, in about 1/2 hour or so.... Looking forward to learning
more about this...