Upper-Triangular Matrix Inverse Using Back-Substitution

[se00an]

Hello all,

Inverting an upper (or lower) triangular matrix is a trivial algorithm, due to the nature of the matrix. I am having an issue getting a part of my upper-triangular matrix inversion function to work, and I would like to get it working soon for a personal project. From reading a verbal description about this special inverse case using back-substitution, I have written some C++ code. I give it an input example 6x6 matrix, it inverts the diagonal correctly, but produces incorrect results above the diagonal. I am checking the true answer using the 'Online Matrix Calculator'. The code is as follows. If you are familiar with this technique, do you have any advice.

If there might be a better forum/thread to post this, I would really appreciate that too, as really want to figure this out.

Thank-you,
Sean Michnowski

// Upper-Triangular Matrix Inversion Test
	double R[36] = {4., 8., 2., 9., 0., 3., 0., 8., 1., 8., 4., 8., 0., 0., 5., 1., 2., 7.,
                    0., 0., 0., 1., 4., 5., 0., 0., 0., 0., 6., 4., 0., 0., 0., 0., 0., 2.};
	double Rinv[36];
	int i, j, k;
	//
	for(i=0; i<36; i++)
		Rinv[i] = 0.0;
	//
	for(j=0; j<6; j++)
	{
		Rinv[6*j+j] = 1./R[6*j+j];
        //
		for(i=0; i<(j-1); i++)
		{
			for(k=0; k<(j-1); k++)
			{
				Rinv[6*i+j] += Rinv[6*i+k]*R[6*k+j];
			}
		}
		//
		for(k=0; k<(j-1); k++)
		{
			Rinv[6*k+j] /= -R[6*j+j];
		}
	}
	//
	displayMatrix("Matrix:", Rinv, 6, 6);
	//
	return 0;

The expected output is:
0.250 -0.250 -0.050 -0.200 0.317 0.667
0.000 0.125 -0.025 -0.975 0.575 0.875
0.000 0.000 0.200 -0.200 0.067 -0.333
0.000 0.000 0.000 1.000 -0.667 -1.167
0.000 0.000 0.000 0.000 0.167 -0.333
0.000 0.000 0.000 0.000 0.000 0.500

Instead, my code produces:
0.25 0 -0.1 -2.25 0.0333333 5.6
0 0.125 0 -1 -0.0833333 2
0 0 0.2 0 -0.0666667 -0.7
0 0 0 1 0 -2.5
0 0 0 0 0.166667 0
0 0 0 0 0 0.5

The diagonals work correctly, but the top is incorrect. Anyone familiar?

[mike_2000_17]

I have to say that your implementation of the backsub looks very bizarre. I have implemented back-substitution is several matrix numerical methods (PLU, QR, HH, Cholesky, SVD, etc.). I have checked them, and they look nothing like yours. Yours looks more like forward-substitution.

First, your bounds at line 14 and 16 are almost certainly wrong. You are definitely skipping elements. I think the bounds of these loops should both be j and not (j-1), because you want to visit the row/col before j, using (j-1) will skip that row/col. And that simply doesn't make sense.

Second, you are, oddly enough, traversing the array in a forward direction while doing back-substitution. You know, it's called back-substitution because you start at the end (last element of the diagonal) and work your way backwards. You are doing the opposite, and I cannot possibly understand how that could still work. If you look, for example, at this site, under the back-sub algorithm, you will notice how the elements are necessarily traversed backwards "(i=n,n-1,n-2,...1)" in the formula.

Finally, the fact that the diagonal is correctly inverted is quite meaningless, the diagonal is almost guaranteed to invert correct if you didn't screw up completely.

PS: Don't ask me to post my back-sub algorithm, it bears no more information than that in the formula at the link given above. And I presume that if you just wanted a ready made algorithm, you would have gotten one already from the numerous places you could get it off the internet.

[se00an]

You are completely right; the j-1 bounds are incorrect. Don't know what I was thinking. That solves it! Thank you so much!