Post Reply 
Problem inverting complex matrix
12-21-2019, 09:48 AM
Post: #1
Problem inverting complex matrix
Consider this complex matrix:

M1 =
\begin{bmatrix}
0.0005 & 0 & 0 & -0.0005 & 0 & 0 & 0 & 1 \\
0 & 9.69696969697E-4+0.00125i & 0 & -0.00125i & -6.66666666667E-4 & -3.0303030303E-4 & 0 & 0 \\
0 & 0 & 3.57142857143E-4-0.0025i & 0 & 0.0025i & 0 & -3.57142857143E-4 & 0 \\
-0.0005 & -0.00125i & 0 & 0.0005+0.00125i & 0 & 0 & 0 & 0 \\
0 & -6.66666666667E-4 & 0.0025i & 0 & 6.66666666667E-4-0.0025i & 0 & 0 & 0 \\
0 & -3.0303030303E-4 & 0 & 0 & 0 & 3.0303030303E-4-1.66666666667E-3i & 0 & 0 \\
0 & 0 & -3.57142857143E-4 & 0 & 0 & 0 & 3.57142857143E-4+0.005i & 0 \\
1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
\end{bmatrix}

Inverting with the HP-Prime using 1/M1 works:
   

Inverting it using inv(M1) fails.
   

Is this a bug or there is something I failed to consider?

If you want to try, M1.hpmat is included in this zip file:

.zip  M1.zip (Size: 289 bytes / Downloads: 8)
Find all posts by this user
Quote this message in a reply
12-21-2019, 12:17 PM
Post: #2
RE: Problem inverting complex matrix
On a 48GX, matrix inversion works ok, and 1/M is not allowed.
Strangely, COND causes a warmstart, but SVD reveals that the matrix is of rank 8.
The smallest singular value is in the order of 2e-4, so the matrix is far from being singular.

Cheers, Werner
Find all posts by this user
Quote this message in a reply
12-21-2019, 01:02 PM
Post: #3
RE: Problem inverting complex matrix
inv calls the CAS, which works better with exact data. You can run approx(inv(exact(...))). If you want me to make a more complete diagnostic, please paste your input as plain text (other formats are more complicated to handle), using the following syntax (replace ... by appropriate data):
m1:=[[0.0005,0,...],[...]]; inv(m1);
Find all posts by this user
Quote this message in a reply
12-21-2019, 01:15 PM
Post: #4
RE: Problem inverting complex matrix
(12-21-2019 01:02 PM)parisse Wrote:  inv calls the CAS, which works better with exact data. You can run approx(inv(exact(...))). If you want me to make a more complete diagnostic, please paste your input as plain text (other formats are more complicated to handle), using the following syntax (replace ... by appropriate data):
m1:=[[0.0005,0,...],[...]]; inv(m1);

Here's his matrix in plain-text format, for easy copy & paste:

M1:=[[0.0005,0,0,-0.0005,0,0,0,1],[0,9.69696969697E-4+0.00125*i,0,-0.00125*i,-6.66666666667E-4,-3.0303030303E-4,0,0],[0,0,3.57142857143E-4-0.0025*i,0,0.0025*i,0,-3.57142857143E-4,0],[-0.0005,-0.00125*i,0,0.0005+0.00125*i,0,0,0,0],[0,-6.66666666667E-4,0.0025*i,0,6.66666666667E-4-0.0025*i,0,0,0],[0,-3.0303030303E-4,0,0,0,3.0303030303E-4-1.66666666667E-3*i,0,0],[0,0,-3.57142857143E-4,0,0,0,3.57142857143E-4+0.005*i,0],[1,0,0,0,0,0,0,0]]

<0|ΙΈ|0>
-Joe-
Visit this user's website Find all posts by this user
Quote this message in a reply
12-21-2019, 04:52 PM
Post: #5
RE: Problem inverting complex matrix
det(exact(m)) returns:
(-1245591211435133830166343338512439683+51639354804136441655581854007061939*i)/44328554791945219310918948796875252340096000000000000, evalf() to -2.80990710679e-17+1.16492303993e-18*i.
The smallest singular values of m is not too small, but there are 6 very small singular values, this is the reason why the determinant is small, and that raises an error at the end of the Gauss algorithm.
As I said, approx(inv(exact(...))) does the job.
Find all posts by this user
Quote this message in a reply
12-21-2019, 08:43 PM
Post: #6
RE: Problem inverting complex matrix
(12-21-2019 04:52 PM)parisse Wrote:  As I said, approx(inv(exact(...))) does the job.

approx(inv(exact(...))) takes over 15 seconds to produce the result.
Find all posts by this user
Quote this message in a reply
12-21-2019, 10:17 PM
Post: #7
RE: Problem inverting complex matrix
Slightly off-topic, but anyone know of a simpler way to start editing a matrix than to type EDITMAT(M1)? I always loved the 50g where a cursor down with matrix on the stack just opened the matrix editor. Even if when a matrix is selected on the screen if there was something like an "EDIT" option on the menu that would launch the matrix editor. As is I find it cumbersome to enter and edit matrices, unless I'm missing something obvious?
Visit this user's website Find all posts by this user
Quote this message in a reply
12-21-2019, 10:42 PM (This post was last modified: 12-21-2019 10:44 PM by jesuscf.)
Post: #8
RE: Problem inverting complex matrix
(12-21-2019 10:17 PM)Jacob Wall Wrote:  Slightly off-topic, but anyone know of a simpler way to start editing a matrix than to type EDITMAT(M1)?

Shift, Matrix, M1.
Shift, Mem, Matrices, M1.
Find all posts by this user
Quote this message in a reply
12-21-2019, 10:57 PM
Post: #9
RE: Problem inverting complex matrix
Quote:Shift, Matrix, M1.

Right, and if I don't want/need a named matrix?

Perhaps even more so for RPN mode, it would be very handy to have:
- Shift Matrix - opens matrix editor
- OK or ENTER pops the matrix on the stack
- From the stack: a) select to reveal EDIT option and b) cursor DOWN to edit

This way to enter two matrices and perform some basic calculations would be IMHO way faster and more intuitive.

Shift Mem Matrices M0-M9 to edit those matrices is fine.
Visit this user's website Find all posts by this user
Quote this message in a reply
12-22-2019, 08:54 AM
Post: #10
RE: Problem inverting complex matrix
(12-21-2019 08:43 PM)jesuscf Wrote:  
(12-21-2019 04:52 PM)parisse Wrote:  As I said, approx(inv(exact(...))) does the job.

approx(inv(exact(...))) takes over 15 seconds to produce the result.

This is expected since computations are done with rationals.
After more tracing with the numeric algorithm, I discovered that it failed because checking of non zero pivot failed: inv(M1) fails while inv(1000*M1) works. This should now be fixed in the source code.
Find all posts by this user
Quote this message in a reply
12-22-2019, 05:35 PM
Post: #11
RE: Problem inverting complex matrix
(12-22-2019 08:54 AM)parisse Wrote:  
(12-21-2019 08:43 PM)jesuscf Wrote:  approx(inv(exact(...))) takes over 15 seconds to produce the result.

This is expected since computations are done with rationals.
After more tracing with the numeric algorithm, I discovered that it failed because checking of non zero pivot failed: inv(M1) fails while inv(1000*M1) works. This should now be fixed in the source code.

Thank you Parisse! I also want to clarify that to get the correct result, inv(1000*M1) needs to be multiplied by 1000: inv(1000*M1)*1000.
Find all posts by this user
Quote this message in a reply
Post Reply 




User(s) browsing this thread: 1 Guest(s)