Precision problem with matrix inversion

[Maro]

when I tried to solve the following original system of linear equations, I observed that obviously the Prime has some precision problems with matrix inversion if the elements have large or small exponents.

\( \underbrace{\begin{bmatrix} 130.2544E9 & 34.4351E9 & 17.1192E9 & 10.4385E9 \\
55.2158E9 & 32.7918E9 & 18.5083E9 & 11.6938E9 \\
26.9576E9 & 29.2367E9 & 22.0630E9 & 15.2900E9 \\
17.8132E9 & 25.8109E9 & 27.0660E9 & 21.6890E9 \end{bmatrix}}_{\textrm{Matrix } P} \cdot
\begin{bmatrix} Q_1 \\ Q_2 \\ Q_3 \\ Q_4 \end{bmatrix} =
\underbrace{\begin{bmatrix} 400E3 \\ 400E3 \\ 400E3 \\ 400E3 \end{bmatrix}}_{\textrm{Vector }V} \)

Solving for vector \([Q_i]\) by simply calculating \([Q_i] = \frac{[V]}{[P]}\) results in

\( \begin{bmatrix} Q_1 \\ Q_2 \\ Q_3 \\ Q_4 \end{bmatrix} =
\begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \end{bmatrix}\)

The correct result should be:

\( \begin{bmatrix} Q_1 \\ Q_2 \\ Q_3 \\ Q_4 \end{bmatrix} =
\begin{bmatrix} -8.611E-9 \\ 6.765E-6 \\ 14.569E-6 \\ -7.782E-6 \end{bmatrix}\)

In contrast to the Prime the 50g gives the correct result without any problem !!

When reducing the exponents on both sides of the matrix equation, for example:

\( \begin{bmatrix} 130.2544E6 & 34.4351E6 & 17.1192E6 & 10.4385E6 \\
55.2158E6 & 32.7918E6 & 18.5083E6 & 11.6938E6 \\
26.9576E6 & 29.2367E6 & 22.0630E6 & 15.2900E6 \\
17.8132E6 & 25.8109E6 & 27.0660E6 & 21.6890E6 \end{bmatrix} \cdot
\begin{bmatrix} Q_1 \\ Q_2 \\ Q_3 \\ Q_4 \end{bmatrix} =
\begin{bmatrix} 400 \\ 400 \\ 400 \\ 400 \end{bmatrix} \)

the result is still a "zero-vector", thus matrix inversion still gives an inverse with 0 for each element.

When solving the system

\( \begin{bmatrix} 130.2544 & 34.4351 & 17.1192 & 10.4385 \\
55.2158 & 32.7918 & 18.5083 & 11.6938 \\
26.9576 & 29.2367 & 22.0630 & 15.2900 \\
17.8132 & 25.8109 & 27.0660 & 21.6890 \end{bmatrix} \cdot
\begin{bmatrix} Q_1 \\ Q_2 \\ Q_3 \\ Q_4 \end{bmatrix} =
\begin{bmatrix} 400 \\ 400 \\ 400 \\ 400 \end{bmatrix} \)

the result from the Prime is correct:

\( \begin{bmatrix} Q_1 \\ Q_2 \\ Q_3 \\ Q_4 \end{bmatrix} =
\begin{bmatrix} -8.611E-3 \\ 6.765 \\ 14.569 \\ -7.782 \end{bmatrix}\)

I tried various settings without any success. When exponents increase the Prime is no longer able to invert the coefficient matrix in a correct way, whereas the 50g always shows the correct results.

Any thoughts about that? Thanks a lot

[roadrunner]

Maro,

I get:

[−8.61069604424e−9,6.76445309046e−6,1.45694987614e−5,−7.78188444368e−6]

Maybe you have a typo?

Road

[roadrunner]

Type it in like this:

   

If you don't type approx you get all zeros.

Road

[Maro]

Thank you Road, you are right. When using approx(...) it works in textbook mode.

But I'm still surprised, since it should work also without using approx(...) in the Home environment, or what do I miss? Is there any reason why the Prime fails without approx(...) in this case? I'm working in RPN mode, putting the two matrices on the stack and divide them ... always resulting in zero-vector. How can I use approx(...) in RPN mode ??

[roadrunner]

I'm sorry Maro, don't have an answer for you. I am not an expert on the prime and I never use it in RPN mode. It's funny, because I never use my HP 50G in algebraic mode.

[Didier Lachieze]

In RPN mode, if M1 is your matrix P and M2 is your vector V, you can do: M1 inv(1) M2 *

[Anders]

It could be that the internal precision is too low when inverting the matrix combined with a multiplication with a vector in RPN mode or simpler: division of a vector with matrix.

"In the good old days...bla bla" for instance with the HP 28s in the 80's, this particular situation was recognized as so common by the HP 28S and HP 48 too I believe (trying to solve a large system of unknowns) so when the user did a division of a vector with matrix from the stack it recognized this situation and used the increased 16 digits precision precisely because this extra precision is needed (if I remember correctly the HP 28s used the standard 12 digits precision for normal calculations).

I do not know if the Prime add extra precision or not for Vector Matrix division, but is should and hopefully this would lead to the right result.

After I wrote this I just had to dig up my HP 28S and guess what - it gets it right!

So I got around to punched in the matrix (let's call it M1) in my Prime too and RPN does not behave as it should with this matrix:
1) V1 M1 / results in the null vector and on top:
2) inverting the matrix: M1 shift / results in the 4x4 null matrix
3) M1 1 / also results in the 4x4 null matrix.
As previously posed
4) M1 inv(1) do result in the correct vector

[Maro]

(07-13-2015 10:40 PM)Didier Lachieze Wrote:  In RPN mode, if M1 is your matrix P and M2 is your vector V, you can do: M1 inv(1) M2 *

thanks Didier, this actually works, but inverting a matrix should also work by simply using the key combination "Shift ÷" (according to the Prime manual!) and in this case it gives the wrong result. Why should we tediously key in "inv(1)" if there is "Shift ÷" on the keyboard.

For me this seems actually to be a bug!

[cyrille de brébisson]

Hello,

I can only speak of the numerical mode here.
The matrix inversion and multiplication have been coded with 3 conflicting constraints:
- High precision
- 0 result elements should be 0
- n^3 (not n^4) and memory constraint

Of course, it is not possible to have all 3 at the same time.
- Calculations ARE done in 15 digit precision on a row by row, but are stored in 12 digit (for memory)
- small elements (smaller in magnitude than epsilon*min number) are transformed into 0
- n^3 algorithm

This leads to the behavior that you observed.
If you do want full precision, that is what the CAS is there for.

Cyrille

[parisse]

I did not check with a CAS numeric matrix, because you did not paste it in algebraic format. As a general rule, if you think you have found a bug, please report it with a textual input that everyone can copy/paste to the emulator, not just as a screen capture. It will cost you more time, but it will save every tester that time, that's much more efficient! Thank you.

[Maro]

(07-15-2015 04:59 AM)cyrille de brébisson Wrote:  Hello,

I can only speak of the numerical mode here.
The matrix inversion and multiplication have been coded with 3 conflicting constraints:
- High precision
- 0 result elements should be 0
- n^3 (not n^4) and memory constraint

Of course, it is not possible to have all 3 at the same time.
- Calculations ARE done in 15 digit precision on a row by row, but are stored in 12 digit (for memory)
- small elements (smaller in magnitude than epsilon*min number) are transformed into 0
- n^3 algorithm

This leads to the behavior that you observed.
If you do want full precision, that is what the CAS is there for.

Cyrille

thanks Cyrille for your comment. Does that mean, that when working in Home environment (RPN):

1) Matrix on the stack: inv(1) ➝ resulting inverse matrix is OK, because inv() is actually a CAS command using full precision
2) Matrix on the stack: Shift ÷ ➝ resulting inverse matrix is NULL-matrix (not OK), because Shift ÷ is using reduced precision of numerical mode

In other words, "inv()" and "Shift ÷" are not the same in Home environment?

[Maro]

(07-15-2015 05:47 AM)parisse Wrote:  I did not check with a CAS numeric matrix, because you did not paste it in algebraic format. As a general rule, if you think you have found a bug, please report it with a textual input that everyone can copy/paste to the emulator, not just as a screen capture. It will cost you more time, but it will save every tester that time, that's much more efficient! Thank you.

Sorry, but I'm not using the emulator as long as there is no proper Mac OSX support. Thus I gave the original matrix in LaTex code (no screen capture!). You should be able to copy the matrix from this code segment ... even without copy/paste it takes about one minute to have the matrix in the calculator ... ?

[Paul Dale]

(07-15-2015 04:59 AM)cyrille de brébisson Wrote:  - n^3 (not n^4) and memory constraint

Just wondering how you get a O(n⁴) matrix inversion algorithm? All the ones I know about are O(n³). Even the more stable ones. Gaussian elimination and the common LU decompositions are O(n³). I believe that there are more efficient algorithms.


- Pauli

[parisse]

(07-15-2015 06:42 AM)Maro Wrote:  Sorry, but I'm not using the emulator as long as there is no proper Mac OSX support. Thus I gave the original matrix in LaTex code (no screen capture!). You should be able to copy the matrix from this code segment ... even without copy/paste it takes about one minute to have the matrix in the calculator ... ?

There is a high risk of error entering data: there are 20 coefficients to enter, it will certainly take at least a couple of minutes. Why not algebraic entry? If someone want to see the matrix, copy/paste it into the emulator or in xcas. LaTeX code looks nice, but it can not be pasted easily.

[Maro]

(07-15-2015 07:46 AM)parisse Wrote:  There is a high risk of error entering data: ...

Oh yes that's absolutely true ... and then you have the basic principle of informatics "shit in ... shit out". Next time I will do it in algebraic mode.

(07-15-2015 07:46 AM)parisse Wrote:  ... there are 20 coefficients ...

16 (4 x 4 matrix)

[parisse]

I said 20 because of the second member :-)

[Anders]

(07-15-2015 04:59 AM)cyrille de brébisson Wrote:  Hello,

I can only speak of the numerical mode here.
The matrix inversion and multiplication have been coded with 3 conflicting constraints:
- High precision
- 0 result elements should be 0
- n^3 (not n^4) and memory constraint

Of course, it is not possible to have all 3 at the same time.
- Calculations ARE done in 15 digit precision on a row by row, but are stored in 12 digit (for memory)
- small elements (smaller in magnitude than epsilon*min number) are transformed into 0
- n^3 algorithm

This leads to the behavior that you observed.
If you do want full precision, that is what the CAS is there for.

Cyrille

So:
1) I suggest this is a bug in the RPN implementation
2) Solving linear equation system through a vector / matrix division is maybe the most basic linear algebraic function that an engineer or a student performs, so it needs to work also in RPN mode
3) HP matrix capable RPN calculators have been able to perform vector - matrix division correctly with sufficiently high precision for ~30 years until now with Prime (I know of HP 28C, HP 28S, 48SX, 50G), so of course it is possible to implement this correctly. These calculators had 2KBs+ of RAM (HP 28C had 2KB) and Prime have 3-4 magnitudes more RAM with 32MB. The "memory constraint" argument does not hold.
4) There are many books and papers written on algorithms to solve vector and matrix division (and matrix inversion) with sufficient precision ranging from O(n^2.3727) to O(n^3). There are also iterative models to improve the result until error < epsilon. In other words, this is a very well research and understood problem. Why would you use an O(n^4) algorithm? (O=Ordo)
5) Suggesting to use CAS mode, which is likely disabled in Exam mode, is not helpful for students if you want them to use RPN (or are you saying we should stop using RPN?)

In general I think HP should thank users (like Maro) for bringing these cases to their attention and address it in an update. You are basically receiving free testing, helping you to make a better product.

[Maro]

(07-15-2015 08:28 PM)Anders Wrote:  ...
2) Solving linear equation system through a vector / matrix division is maybe the most basic linear algebraic function that an engineer or a student performs, so it needs to work also in RPN mode
....

This is basically the point. And it needs to work properly in algebraic and textbook mode too, which is NOT the case at the moment. It only works fine when using "approx([V]/[P])" but not for simply doing "[V]/[P]", which would be the most obvious way when working in Home environment ... Why that??

(07-15-2015 08:28 PM)Anders Wrote:  ...
3) HP matrix capable RPN calculators have been able to perform vector - matrix division correctly with sufficiently high precision for ~30 years until now with Prime (I know of HP 28C, HP 28S, 48SX, 50G), so of course it is possible to implement this correctly. These calculators had 2KBs+ of RAM (HP 28C had 2KB) and Prime have 3-4 magnitudes more RAM with 32MB. The "memory constraint" argument do not hold.
...

This is another point. The "old fashioned" 50g can do it without any problem in numerical mode (exact mode is not needed!). So what is the reason why the "powerful" Prime shouldn't be able to do this matrix inversion in a correct way in numerical mode too?? Perhaps somebody involved in the Prime development can further explain why the 50g succeeds and the Prime fails? What is the difference?

(07-15-2015 08:28 PM)Anders Wrote:  ...
5) Suggesting to use CAS mode, which is likely disabled in Exam mode, is not helpful for students if you want them to use RPN (or are you saying we should stop using RPN?)
...

+1

[Thomas Radtke]

(07-15-2015 08:28 PM)Anders Wrote:  5) Suggesting to use CAS mode, which is likely disabled in Exam mode, is not helpful for students [...]

I hope exams are about e.g. doing matrix inversions by hand, not proving knowledge in using a calculator ;-).

[Werner]

Cyrille:

The 12/15 digit algorithm as used by the venerable 42S has no problems with this system of equations, so I don't quite see why the Prime would.

Everyone else:
B A / is not the same as inv(A)*B.. The former is faster and more accurate, as it does not compute the inverse. Or at least that's how it has been in the 71B, 42S, 48/49/50.

Cheers, Werner