Calculating an inverse matrix in Matlab

Question

I'm running an optimization algorithm that requires calculation of the inverse of a matrix. The goal of the algorithm is to eliminate negative values from the matrix A and obtain the new matrix B. Basically, I start with known square matrices B and C of the same size.

I start by calculating the matrix A which is equal to:

A = B^-1 * C

Or in Matlab:

A = B\C;

I use this because Matlab told me B\C is more accurate than inv(B)*C.

The negative values in A are then divided by two and A is then normalised so that it's rows have length of 1. Using this new A, I calculate a new B with:

(1/N) * A * C' = B^-1

where N is just a scaling factor (# of columns in A). This new B would then be used again in the first step and these iterations continue until the negatives in A are gone.

My problem is I have to calculate B from the second equation and then normalise it.

invB = (1/N)*A*C';
B = inv(invB);

I've been calculating B using inv(B^-1) but after a few iterations I start getting messages that B^-1 is "close to singular or badly scaled."

This algorithm actually works for smaller matrices (around 70x70) but when it gets up to about 500x500 I start getting these messages.

Are there any better ways to calculate inv(B^-1)?

Answer 1

You should definitely head warnings about singular matrices. Results in numerical linear algebra tend to break down as you move toward matrices with high condition numbers. The underlying idea is if

A*b_1 = c

and we're actually solving the problem (because we are using approximate numbers when we use computers)

(A + matrix error)*b_2 = (c + vector error)

how close are b_1 and b_2 as a function of the matrix and vector errors? When A has small condition number b_1 and b_2 are close. When A has large condition number b_1 and b_2 are not close.

There is an informative piece of analysis you could do on your algorithm. At each iteration, after you've found B, find use Matlab to find the condition number of it. This is

cond(B)

You will likely see the number climb rapidly. This indicates that every time you iterate your algorithm, you should trust your result for B less and less.

Problems like this crop up all the time in numerical mathematics. If you'll be working with numerical algorithms frequently you should take some time to familiarize with the role of condition numbers in the field and preconditioning techniques as mentioned above. My preferred text for this is "Numerical Linear Algebra" by Lloyd Trefethen, but any text on Numerical Algebra should address some of these issues.

Best of luck, Andrew

Answer 2

The main issue is that your matrix has a high condition number (i.e. really small rcond(B) in your case). This is due to the iterative structure in your algorithm, I guess. As you do each iteration your small singular values get smaller and smaller so your condition number grows exponentially. You should check preconditioning to avoid this kind of behavior.