In the previous lab, we introduced the idea of solving a single linear system A*x=b by constructing the PLU factorization of A.
We will now look at the PLU factorization more closely. In particular, we will see that we can use this tool to solve multiple linear systems, to compute the determinant or inverse, or to solve linear systems involving the transposed matrix. These applications make the PLU factorization a very useful tool.
The PLU factorization is a relatively expensive operation, and it would be worth our while to consider whether we can use the factors for other purposes. In fact, one of the simplest reuses involves the solution of multiple linear systems, where the system matrix is fixed, but we have several right hand sides.
If we've just solved A*x1=b1, then we actually did two things:
[ P, L, U ] = ge_pp ( A )
x1 = plu_solve ( P, L, U, b1 )
x2 = plu_solve ( P, L, U, b2 )
Exercise - use one factorization to solve two systems. Use the Frank matrix of size n=5, and set
A = frank(n)Factor A once, and solve linear system 1, and then solve linear system 2. I was sneaky in defining b2, but if I point out that I'm setting it to column n-2 of A, can you figure out, beforehand, what the solution x2 must be?
x1 = [1 : n]'
b1 = A * x1
b2 = A(:,n-2)
If both right hand sides are available at one time, MATLAB would actually prefer to solve them together. To do this, we would have to set up a single N by 2 array for the right hand side, and make sure our plu_solve routine was written properly to handle simultaneous multiple right hand sides.
M File - modify your plu_sol routine so that it can handle multiple right hand sides, stored in a single array b. Here are some suggestions:
[ n, nrhs ] = size ( b ) ... (figure out number of RHS's)
z = P' * b ... (can stay the same!)
y = zeros ( n, nrhs )
for k = 1 : nrhs
for j = 1 : n
y(j,k) = ?
for i = j+1 : n
z(i,k) = z(i,k) - L(i,j) * y(j,k);
end
end
end
...similar changes for the U code.
Exercise - test your modified solver. Use the dif2 matrix of size n=5, and set
A = dif2(n)Factor A once, and then solve both linear systems by a single call to plu_sol.
b(:,1) = A * [1:n]'
b(:,2) = A(:,n-2)
The determinant is easy to compute because of the following facts:
det ( A ) = det ( P ) * det ( L ) * det ( U )
= (+1 or -1) * U11 * U22 ... * Unn.
M File - Write a routine called plu_det.m which computes the determinant of a matrix given its PLU factorization. The routine should have the form:
det = plu_det ( P, L, U )Copy the M file p_det.m from the web page, which will compute the determinant of P. You can either call p_det from within your routine, or simply copy the code into your routine. You should be able to write the rest of the routine yourself.
Exercise - Use your routine plu_det.m to compute the determinants of a few matrices:
Matrix N Determinant plu_det(A)
DIF2 5 -6.0 ___________
DIF2 20 21.0 ___________
Frank 5 1.0 ___________
Hilbert 5 3.7E-12 ___________
We could try to compute the inverse of A by computing the inverses of each of the PLU factors and multiplying them out. (This isn't too hard.) But a better way, which is faster, and re-uses code we've already written and tested, is to construct the inverse by solving the following set of linear systems:
A * X = Iwhere I is the identity matrix. (The identity matrix of order n is called eye(n) in MATLAB.) The solution matrix X will actually be the inverse. Because we've fixed up our plu_solve routine, we can do this in a single operation.
Exercise - Create a routine called inverse.m which computes the inverse of a matrix. The routine should have the form
function B = inverse ( A )Inside of your routine, what do you have to do in order to compute the inverse? Test your routine on the 5 by 5 dif2 matrix, and compare your results to the output of the official MATLAB inv command, or the dif2_inv routine I gave you.
I said that if you had the PLU factors of a matrix A, you could also use them to solve linear systems involving the transpose matrix, that is, equations of the form:
A' * x = bIn the assignment, I'm going to ask you to try to write the code to do this, so pay attention!
Consider the fact that if
A = P * L * Uthen the transpose matrix has the factorization
A' = U' * L' * P'If we needed to solve systems involving A', we could go ahead and compute the PLU factorization of A', but with a little work, we can use this U'L'P' factorization instead.
(PLU)' Solution Algorithm: To solve the transposed linear system A'*x=b, given factors P, L, U (which are factors for the untransposed matrix) and right hand side b,
Assignment - Make a copy of your code plu_solve.m, calling it plut_solve.m. Modify this code so that it can solve a linear system involving the transposed matrix, given PLU factors of the original matrix.
Design suggestions: - To start this code, simply take the three sections of your original program, and reverse their order. Now rewrite each section so that it's solving the transposed equation. Start with the transposed P system, which is easy. To solve the transposed L system, you have to realize that L' is a unit upper triangular matrix, and that its (I,J) entry is stored in L(J,I). So essentially, rewrite the code to look like your upper triangular solve code, but adjust for the fact the your matrix now has a unit diagonal, and the actual elements have the indices reversed. Finally, try to set up the U' solver, which looks like your old L solver, but now the lower triangular matrix doesn't have a unit diagonal, and, again, the entry indices are flipped.
Test: When you think your code is working, carry out the following test:
n = 5(Do not try your code on the Hilbert or dif2 matrices. They are bad examples for this problem, because they are symmetric!)
A = frank ( n )
[P,L,U] = ge_pp ( A )
x = [1:n]'
b = A' * x
x2 = plut_solve ( P, L, U, b )
Results: If your code has passed the test I just descrbed, it is probably working correctly. I want to take a look at your code (that's all). Mail your plut_solve.m code to me.