LU Decomposition Solution

In this assignment, you will implement a Matlab function to decompose a matrix into lower and upper triangular matrices (L and U), i.e., PA = LU where P is a row permutation matrix, and apply it to solve a computational physics problem.

1 Download


For Section 6, we provide codes that can compute forces in a simple physical system and visualize the results. This code has recently been updated. Download the file sec6code.zip from the Moodle submission page, and extract it into the folder where you will be working on the assignment. You should obtain the files get_forces.m, draw.m, and a directory +defo.

2 Submission Guidelines


You will submit a zip file that contains the following .m files on Moodle:

• replacement_lu.m

• interchange_lup.m

• my_lu.m

• my_lup.m

• build_matrix.m

• get_displacements.m

Note 1: Each file must be named as specified above. Note 2: Each function must return the specified value.
Note 3: You may use Matlab’s high-level array manipulation syntax, e.g. size(A), A(i,:), and [A,B], but the built-in linear algebra functions such as inv, lu, rref, and A\b are not allowed (except in Section 6). Please contact the TAs for further questions.

Note 4: No collaboration is allowed. Plagiarism cases will be directly reported to the University.





3 Elementary Row Operations (4 points)


You will implement two simple modular functions that will be useful for LU decomposition. These are very similar to the functions you implemented in Assignment 1.

Notation: for subsequent sections, Ai indicates the ith row of A, and Aij indicates the (i, j) entry of A.

Specification:

function [U_new, L_new] = replacement_lu(U, i, j, s, L)
Input: two square matrices U and L, two integers i and j, and a scalar s.
Output:

• Unew : updated U by subtracting s times row j from row i, i.e. performing the row replacement operation Ui ← Ui − sUj .

• Lnew : updated L by filling in its (i, j) entry with s, i.e., Lij ← s.

Algorithm 1 replacement lu
1: n ← the number of columns of U
2: Lij ← s
3: for 1 ≤ k ≤ n do
4: Uik ← Uik − sUjk . Row replacement
5: end for
6: return L, U

Warning! If you are adapting your code from Assignment 1, be careful that we are now subtracting s times row j from row i instead of adding it.





function [U_new, L_new, P_new] = interchange_lup(U, i, j, L, P) Input: three same size square matrices U, L, and P, and two integers i and j. Output:

• Unew : updated U by swapping rows i and j, i.e. Ui ↔ Uj .

• Lnew : updated L by swapping only the first (i − 1) entries of rows i and j, i.e., Li,1 ↔ Lj,1 , . . . , Li,(i−1) ↔ Lj,(i−1) as shown in Figure 1.

• Pnew : updated P by swapping rows i and j, i.e. Pi ↔ Pj .

i-1 i

i-1 i
1 0 0  0 

1 0 0  0 
∗  

∗  
   
i ∗ ∗ 1
∗ ∗ 0 

Partial row
  exchange

i ∗ ∗ 1  
∗ ∗ 0  
   
j ∗ ∗
∗ ∗

0  1 0 

0  
n

j ∗ ∗
∗ ∗

0  1 0 

0  

Figure 1: Partial row exchange in L.


Algorithm 2 interchange lup
1: Ui ↔ Uj . Row interchange
2: Pi ↔ Pj . Row interchange
3: if i 1 then
4: for 1 ≤ k ≤ i − 1 do
5: Lik ↔ Ljk . Partial row interchange
6: end for
7: end if
8: return L, U , P





4 LU Decomposition (4 points)


In this part, you will write a function that performs LU decomposition without pivoting. We will deal with pivoting in the next part of the assignment.

Specification:

function [L, U] = my_lu(A) Input: an n × n square matrix A. Output:

• L: an n × n lower triangular matrix where the diagonal entries are all one,
 1 0 0 
e.g.,  ∗ 1 0  where ∗ is a potentially nonzero element.
∗ ∗ 1

• U: an n × n upper triangular matrix.

To get full credit, your function should handle the following case:

• Early termination: Due to round-off error, there may exist free variables whose pivots are extremely small but not precisely zero. You should terminate your LU decomposition if the absolute value of a pivot is less than 10−12 .

The process of LU decomposition uses Gaussian elimination that transforms A to an upper triangular matrix U while recording the pivot multipliers in a lower triangular matrix L.

1. Initialize L to the identity matrix, and U to A. You can use Matlab’s built-in function eye(n).
2. At the ith step, for each row k below the ith row,

(a) Record the pivot multiplier to L at (k, i), i.e., Lk,i = Uk,i /Ui,i as shown in Figure 2. Note that this fills in the ith column of L.

(b) Reduce the kth row of U using the pivot multiplier, i.e., Uk =
Uk − (Uk,i /Ui,i )Ui where Ui is the ith row.




∗ 

 0  
   
i ∗ ∗ 1

  U = i

0 0

Uii  
L =    
∗ ∗   

0 0 0  
k ∗ ∗

U / U

 1 0 

k 0 0

0  ∗ ∗ 
 ki ii   
∗ ∗

  

  

  


Figure 2: LU decomposition at the ith step.





We provide pseudocode for LU decomposition in Algorithm 3.

Algorithm 3 LU decomposition of A
1: n ← the number of columns of A
2: L ← In where In is n × n identity matrix.
3: U ← A
4: for 1 ≤ i ≤ n − 1 do
5: if Uii is nearly zero* then
6: return L, U . Early termination
7: end if
8: for i + 1 ≤ k ≤ n do
9: p ← Uki /Uii
10: [U , L] = replacement_lu(U , k, i, p, L)
11: end for
12: end for
13: return L, U

*Consider a number to be nearly zero if its absolute value is smaller than 10−12 .

Note: When terminating early, L is not unique, i.e., the order of rows correspond- ing to zero rows in U can be different. As long as the resulting decomposition satisfies A = LU, the solution L and U are valid.

Test cases:

• [L,U] = my_lu([4,-2,2;-2,5,3;2,3,9])
 1 0 0

4 −2 2
L = −0.5 1 0, U = 0 4 4
0.5 1 1

• Early termination:

0 0 4
[L,U] = my_lu([4,-2,2,1;-2,5,3,3;4,-2,2,1;-2,5,3,3])
 1 0 0 0

4 −2 2 1 
L = −0.5 1 0 0, U = 0 4 4 3.5.
 1 0 1 0

0 0 0 0 
 
−0.5 1 0 1

 
0 0 0 0
Without early termination, it will return L and U with NaN (Not-a-Number)
entries.

• You can test your algorithm by constructing a random n × n matrix A, e.g., A=rand(10,10), and testing whether multiplying L and U reconstructs A. You also need to check whether L is a lower triangular matrix with diagonal entries equal to 1, and U is an upper triangular matrix.





5 LU Decomposition with Partial Pivoting
(4 points)


Based on your my_lu, you will write numerically stable LU decomposition with partial pivoting. At the ith step of LU decomposition (ith pivot column), you will find the row that has the largest absolute value in the pivot column (say row j), and swap the ith and jth rows of U as usual. Simultaneously, you will swap the partial entries of the ith and jth rows of L, and record the row exchange in a permutation matrix P. For further details, please see http://www.math.kent.edu/~reichel/courses/intr.num.comp.1/fall09/lecture9/lecture4.pdf

Specification:

function [L, U, P] = my_lup(A) Input: an n × n square matrix A. Output:

• L: an n × n lower triangular matrix where the diagonal entries are all 1.

• U: an n × n upper triangular matrix.

• P: an n × n permutation matrix.

The process of LU decomposition with partial pivoting needs to compute an additional row permutation matrix P.

1. Initialize L and P to the identity matrix, and U to A. You can use
Matlab’s built-in function eye(n).

2. At the ith step,

(a) Similar to Assignment 1, perform partial pivoting in U.

(b) Record the row exchange to the permutation matrix P by exchanging its corresponding rows.

(c) Exchange the corresponding row entries in L to reflect the row ex- change in U. Note that this exchange only involves with the row entries that have been recorded, i.e., not entire row exchange.

(d) For each row k below the ith row, record the pivot multiplier to L and replace the row in U using the pivot multiplier, like in the previous part of this assignment.





We provide pseudocode for LUP decomposition in Algorithm 4.

Algorithm 4 LUP decomposition of A
1: n ← the number of columns of A
2: L ← In where In is n × n identity matrix.
3: P ← In . Permutation initialization
4: U ← A
5: for 1 ≤ i ≤ n − 1 do
6: j ← the row index of the largest absolute value among rows i to n in the
ith column of U
7: [U , L, P ] = interchange_lup(U , i, j, L, P )
8: if Uii is nearly zero* then
9: return L, U , and P . Early termination
10: end if
11: for i + 1 ≤ k ≤ n do
12: p ← Uki /Uii
13: [U , L] = replacement_lu(U , k, i, p, L)
14: end for
15: end for
16: return L, U , P


The red lines specify the major modification for partial pivoting.
*As before, consider a number to be nearly zero if its absolute value is smaller than 10−12 .

Note: When terminating early, L and P are not unique, i.e., the order of rows corresponding to zero rows in U can be different. As long as the resulting decomposition satisfies PA = LU, the solution L and P are valid.

Test cases:

• Partial pivoting:
[L,U,P] = my_lup([-2,5,3;2,3,9;4,-2,2])
 1 0 0

4 −2 2 

0 0 1
L =  0.5 1 0, U = 0 4 8 , and P = 0 1 0
−0.5 1 1

• Early termination:

0 0 −4

1 0 0
[L,U,P] = my_lup([4,-2,2;-2,5,3;4+5E-13,-2+2E-13,2+7E-13])
 1 0 0

4 −2 2

0 0 1
L = −0.5 1 0, U = 0 4 4, and P = 0 1 0
1 0 1

0 0 0

1 0 0

• You can test your algorithm by comparing its results with Matlab’s built-in function, [l, u, p] = lu(A), on a randomly generated n × n matrix A, e.g., A=rand(10,10).





6 Force-Displacement Computations (3 points)


Numerical linear algebra is used extensively in computational physics, with applications in a wide range of fields including science, engineering, robotics, and animation. Often, we represent a physical object as a collection of vertices, and we have a computational model that can predict the internal forces that would result from any given deformation of the object (i.e. any given displacement of the vertices). Linear algebra allows us to do the reverse: apply any specified external forces on the object, and predict how it will deform as a result.


1.5

1.5

1 1

0.5

0.5

0 0

-0.5


-0.5 0 0.5 1 1.5 2 2.5 3 3.5

-0.5



-0.5 0 0.5 1 1.5 2 2.5 3 3.5


Figure 3: Left: The rest shape of a simple system with four free vertices. Right: A displacement of the vertices, and the resulting internal forces.

For example, consider the object above which has four free vertices (in black) that can move around, and two constrained vertices (in red) that remain fixed. We can collect the 2D displacements of the 4 free vertices into a single 8-dimensional vector d. We have provided for you a function f = get_forces(d) that takes this vector of displacements d as input and returns the resulting forces f on the vertices, similarly packed into an 8-dimensional vector. In fact, the forces are linear in the displacements, so the function is really a linear transformation get forces : R8 → R8 , and there exists some 8 × 8 matrix A such that f = Ad.

To visualize this relationship, we have provided a function draw(), which can ei- ther be called without arguments to draw the rest shape, or with two arguments draw(d,f) to visualize a deformed shape with forces as arrows. Try choos- ing a vector d ∈ R8 with small entries and calling draw(d, get_forces(d)) to produce something like Fig. 3 (right). (Actually, you may have to draw
0.1*get_forces(d) because the forces are very large.)





Your task is to implement a function d = get_displacements(f, ...) that computes the inverse of this transformation: given the internal forces on the free vertices, find the displacements that would give rise to them. To do this, you will have to find the matrix A for the transformation, compute its LUP factorization, and then use the factorization to quickly solve Ad = f to get d for any given f .

Specification:

A = build_matrix()
Input: None.
Output: the 8 × 8 matrix A corresponding to the linear transformation get_forces. For any vector d ∈ R8 , the matrix-vector product A*d should equal get_force(d).
Hint: Construct the vector e1 = [1, 0, . . . , 0]T and call get_forces on it. What you get must be the first column of A. Do the same for the rest of the columns.

Once you have A, use either my_lup or Matlab’s lu to compute its LUP
factorization. You will need it for the following function:

d = get_displacements(f, L, U, P)
Input: a vector of forces f ∈ R8 , and the LUP decomposition of A.
Output: the vector of vertex displacements d that result in these internal forces, i.e. the solution of Ad = f .

Since PA = LU, we can write y g
PAd = Pf ⇐⇒ L zU}|d{ = zP}|f{ .

Thus we can compute d in three steps: compute the vector g = Pf , solve Ly = g for y, and then solve Ud = y for d. In this function, you may use Matlab’s backslash operator A\b to perform the two solves, since Matlab automatically performs backsubstitution if the matrix is triangular.

Note: In the original code we provided, it turned out that the system matrix you obtained from build_matrix did not require any pivoting, so you did not actually need P to get the right answer. We have now provided an updated version of the code which does induce pivoting. (It still models the same system, but the ordering of the entries in f is different.) Download the updated code and verify that your solution still works on it.


Once you have these functions working, you can compute the deformation caused by any external forces fext : The resulting displacement is the one where the inter- nal and external forces cancel out, so d = get_displacements(-f_ext, ...) with a negative sign. For example, choose fext = [0, 0, 0, 0−0.2, −0.2, −0.2, −0.2]T to apply a downward gravity force to all vertices, compute d, and visualize the result using draw(d, f_ext).





Test cases:

These have been updated for the new version of the provided code.

For verification, we give the top left 2 × 2 block of the matrix A you should obtain:
−23.5355 −3.5355 · · ·
A ≈ 

0 0 · · ·
 . .

. . .

Pulling the fourth vertex upwards with 1 unit of force, that is, f_ext = [0;0;0;0; 0;0;0;1], should produce a displacement d = get_displacements(-f_ext, L,U,P) of approximately
 0.0185 
 0.1432 
 
−0.0269
 0.1212 
d ≈ 0.0149 .
 
 
 0.1950 
 
 0.0177 
 
0.2063