# Diagonal of inverted matrix

95 views (last 30 days)
Hildo on 18 Mar 2017
Edited: James Tursa on 21 Mar 2017
I am working with a large sparse matrix. How can I get the diagonal of the inverse of this matrix? If I use diag(inv(A)) , returns some warnings. Is this the better way?

John D'Errico on 18 Mar 2017
Is that the right way to do it? Very possibly there is no good way. If you are getting warning messages, that generally means your matrix is numerically singular. Finding the inverse of a numerically singular matrix will not be well posed, no matter what computation you use.
So the very first thing you need to do is test the condition number of the matrix. If it is truly very large and sparse, then condest may be the best tool, to give at least an estimate of the approximate condition number. Of course, a lot of people think their matrices are large and sparse, when they are neither truly large or truly sparse. So cond may suffice for you, to tell you if the matrix is singular. Again, if your matrix is singular, then you are wasting your time to compute the diagonal of the inverse, since the inverse matrix will be numerical garbage.
A better solution may depend on how the matrix was created, using a little mathematics. But that is something we are not able to know, since you have told us nothing of value.

Hildo on 20 Mar 2017
The origin of the matrix if some admittance matrix of a electrical system (we call Ybus). A have to get Zbus (Zbus=Ybus^-1). But the important to me is just the diagonal of Zbus. Annex same sample data.
John D'Errico on 20 Mar 2017
As is often the case, people think they have large sparse matrices, when they don't. :)
No matter how sparse it is, a 22x22 matrix is not large. Not even worth using sparse storage to store it. The point is, just make it a full matrix. Things get easier then.
YB = full(Ybus_pu);
>> rank(YB)
ans =
21
>> cond(YB)
ans =
3.22109878812164e+17
The matrix is singular. Many people don't understand what that means. The condition number is roughly 3e17. What that means is if you try to solve a linear system of equations, OR compute the inverse matrix, the system will amplify any noise in your problem by roughly a factor of 3e17.
Is there noise in your problem? YES, there is! The noise comes from how those numbers are stored. They have random junk in the least significant bits of the numbers. That last bit will be corrupted, even if the numbers themselves were computed with no "error". And ANY floating point computations end up corrupting those least significant bits. Even just the process of solving for the inverse.
So accept that there is junk in your matrix entries down in the least significant bits, that is on the order of eps*YB(i,j). In double precision, eps is:
eps
ans =
2.22044604925031e-16
But remember that the condition number of your matrix, thus the extent of any amplification of the noise, is 3e17.
What does this tell you? It says that the elements of the inverse are complete junk. They will be completely corrupted by the noise in those least significant bits of the matrix.
I know. you don't believe me. Lets do a little test. First, compute the diagonal elements of the inverse matrix directly.
diag(inv(YB))
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.161271e-18.
ans =
-6991140573317.74 + 2649900922810.41i
-6991140537171.36 + 2649900963812.15i
-6991140537171.34 + 2649900963812.15i
-6991140537171.26 + 2649900963812.17i
-6991140537171.32 + 2649900963812.16i
-6991140537171.28 + 2649900963812.18i
-6991140537171.24 + 2649900963812.2i
-6991140537171.21 + 2649900963812.21i
-6991140537171.22 + 2649900963812.21i
-6991140537171.18 + 2649900963812.22i
-6991140537171.3 + 2649900963812.17i
-6991140537171.29 + 2649900963812.18i
-6991140537171.21 + 2649900963812.2i
-6991140537171.27 + 2649900963812.19i
-6991140537171.25 + 2649900963812.2i
-6991140537171.24 + 2649900963812.22i
-6991140537171.21 + 2649900963812.23i
-6991140537171.24 + 2649900963812.2i
-6991140537171.23 + 2649900963812.23i
-6991140537171.21 + 2649900963812.24i
-6991140537171.18 + 2649900963812.25i
-6991140573317.74 + 2649900922810.41i
Now, perturb the matrix elements by a TINY amount, on the order of eps.
diag(inv(YB + randn(size(YB)).*eps(YB)))
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.469158e-18.
ans =
1341274048446.89 + 46993008610.608i
1341274044976.97 + 46993074953.964i
1341274044976.98 + 46993074953.9822i
1341274044977.06 + 46993074953.9981i
1341274044976.99 + 46993074954.0002i
1341274044977.03 + 46993074954.0155i
1341274044977.06 + 46993074954.0322i
1341274044977.09 + 46993074954.0497i
1341274044977.09 + 46993074954.0424i
1341274044977.12 + 46993074954.0529i
1341274044977.01 + 46993074954.0161i
1341274044977.02 + 46993074954.0311i
1341274044977.09 + 46993074954.0471i
1341274044977.03 + 46993074954.0459i
1341274044977.04 + 46993074954.0592i
1341274044977.04 + 46993074954.0724i
1341274044977.08 + 46993074954.0829i
1341274044977.06 + 46993074954.057i
1341274044977.06 + 46993074954.0869i
1341274044977.08 + 46993074954.0947i
1341274044977.1 + 46993074954.1026i
1341274048446.89 + 46993008610.608i
The tiny permutations in those elements results in crap that was as large as the elements of the original inverse.
Essentially, if you think of this as a signal to noise thing, there is NO signal remaining in the elements of that inverse.
It does not matter how much you want to compute the elements of the inverse matrix when it is singular. The numbers you will produce are COMPLETELY MEANINGLESS. There is no information content remaining.
Ok, some might now say, but a 22x22 matrix is small. Just use the symbolic toolbox. The condition number is still 3e17. And the noise in your elements is of the same magnitude, because they are created in double precision. That means you will still see amplification of that noise by roughly the condition number.
Sorry, but you can't succeed via that route.
Ok, suppose you go back to the original matrix, and created it in full symbolic form. So never go through double precision. Now all the entries of the matrix are symbolic, and have no corruption in the least significant bits. Can we possibly now survive?
The question is why is your matrix singular. If I look at the singular values, of this thing, I see what is one effectively zero singular value. (Compare it to the largest singular value. It is relatively near eps.)
svd(YB)
ans =
5235.34878600224
224.456140741978
208.213133465402
183.772675608031
144.79532279094
122.428521425073
117.857092301615
100.073512943469
71.6504865943657
66.4367745517612
49.3938713012852
43.2493963956952
35.2213230481699
21.088537030153
18.9310735629174
16.497370799831
14.506749928581
8.95375281077995
4.90625096086107
2.01074591372002
9.2676852076552e-06
1.62533009087101e-14
So I have no idea how that matrix was generated. It may well be that even if you built it in symbolic form, it would still be singular!

Walter Roberson on 18 Mar 2017

Walter Roberson on 20 Mar 2017
If you have the symbolic toolbox, then you can proceed symbolically:
dinv = diag( inv( sym(Ybus_pu) ) );
double(dinv)
The values are mostly close to -8327187525072.06 + 2366252476427.26i with the "ones" and the decimals varying -- the first 12 places are pretty constant for most of the entries.
Walter Roberson on 20 Mar 2017
If you take
Y = sym(Ybus_pu);
Y1 = Y;
Y1(1,1) = Y1(1,1) + 8.11130830789689e-14;
dinv1 = diag(inv(Y1));
r41 = real(dinv1(4));
Y2 = Y;
Y2(1,1) = Y2(1,1) + 1.41747416292681e-13;
dinv2 = diag(inv(Y2));
r42 = real(dinv2(4));
then r41 will be about -15814803937051 and r42 will be about 15828133351471 . This indicates that a change of 1E-14 to 1E-13 can change the sign of the result completely.
But
eps(real(Ybus_pu))
is 4.54747350886464e-13 . which is about 4 to 8 times larger than those shifts.
This tells us that the answers you get out through the process are essentially numeric garbage, completely different with a variation in values in the input smaller than MATLAB double precision can represent.
Your situation is hopeless unless you can generate those bus values to higher precision such as by using the Symbolic Toolbox when you create them.