Diagonal of inverted matrix

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?

Answer 1

[John D'Errico]

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]

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]

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!

Answer 2

[Walter Roberson]

See the discussion at http://www-users.cs.umn.edu/~saad/PDF/Sparse_days_06_15_2010.pdf

[Walter Roberson]

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]

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.