Solving Linear System of Equations The Undo button for Linear - - PowerPoint PPT Presentation

Solving Linear System of Equations

The “Undo” button for Linear Operations

Matrix-vector multiplication: given the data 𝒚 and the operator 𝑩, we can find 𝒛 such that 𝒛 = 𝑩 𝒚 What if we know 𝒛 but not 𝒚? How can we “undo” the transformation?

𝒚 𝒛 𝑩

transformation

𝑩!𝟐 𝒛 𝒚 ?

Solve 𝑩 𝒚 = 𝒛 for 𝒚

Image Blurring Example

• Image is stored as a 2D array of real numbers between 0 and 1

(0 represents a white pixel, 1 represents a black pixel)

• 𝒚𝒏𝒃𝒖 has 40 rows of pixels and 100 columns of pixels
• Flatten the 2D array as a 1D array
• 𝒚 contains the 1D data with dimension 4000,
• Apply blurring operation to data 𝒚, i.e.

𝒛 = 𝑩 𝒚 where 𝑩 is the blur operator and 𝒛 is the blurred image

Blur operator

𝒛 = 𝑩 𝒚

"original” image (4000,) blurred image (4000,) Blur operator (4000,4000) Blur operator

𝒚 𝒛

𝑩

”Undo” Blur to recover original image

Solve 𝑩 𝒚 = 𝒛 for 𝒚 Assumptions:

• 1. we know the blur
• perator 𝑩
• 2. the data set 𝒛 does not

have any noise (“clean data” 𝒚 𝒛

”Undo” Blur to recover original image

Solve 𝑩 𝒚 = 𝒛 for 𝒚

𝒛 + 𝑏 ∗ 10!" (𝑏 ∈ 0,1 ) How much noise can we add and still be able to recover meaningful information from the original image? At which point this inverse transformation fails? We will talk about sensitivity of the “undo” operation later. 𝒛 + 𝑏 ∗ 10!# (𝑏 ∈ 0,1 )

Linear System of Equations

We can start with an “easier” system of equations… How do we actually solve 𝑩 𝒚 = 𝒄 ? Let’s consider triangular matrices (lower and upper):

𝑀!! 𝑀"! 𝑀"" … … ⋮ ⋮ 𝑀#! 𝑀#" ⋱ ⋮ … 𝑀## 𝑦! 𝑦" ⋮ 𝑦# = 𝑐 𝑐" ⋮ 𝑐# 𝑉!! 𝑉!" 𝑉"" … 𝑉!# … 𝑉"# ⋮ ⋮ ⋱ ⋮ … 𝑉## 𝑦! 𝑦" ⋮ 𝑦# = 𝑐 𝑐" ⋮ 𝑐#

2 3 2 1 2 1 3 6 4 2 𝑦! 𝑦" 𝑦# 𝑦$ = 2 2 6 4

Example: Forward-substitution for lower triangular systems

2 𝑦$ = 2 → 𝑦$= 1 3 𝑦$ + 2 𝑦% = 2 → 𝑦%= 2 − 3 2 = −0.5 1 𝑦$ + 2 𝑦% + 6 𝑦& = 6 → 𝑦&= 6 − 1 + 1 6 = 1.0 1 𝑦$ + 3 𝑦% + 4 𝑦& + 2 𝑦# = 4 → 𝑦&= 4 − 1 + 1.5 − 4 2 = 0.25

𝑦! 𝑦" 𝑦# 𝑦$ = 1 −0.5 1.0 0.25

Triangular Matrices

𝑉!! 𝑉!" 𝑉"" … 𝑉!# … 𝑉"# ⋮ ⋮ ⋱ ⋮ … 𝑉## 𝑦! 𝑦" ⋮ 𝑦# = 𝑐 𝑐" ⋮ 𝑐#

𝑦! 𝐕 : , 1 + 𝑦" 𝐕 : , 2 + ⋯ + 𝑦% 𝐕 : , 𝑜 = 𝒄 Hence we can write the solution as

𝑉'' 𝑦' = 𝑐'

Recall that we can also write 𝑽 𝒚 = 𝒄 as a linear combination of the columns of 𝑽

𝑦$ 𝐕 : , 1 + ⋯ + 𝑦'!$ 𝐕 : , 𝑜 − 1 = 𝒄 −𝑦' 𝐕 : , 𝑜 → 𝑉'!$,'!$ 𝑦'!$ = 𝑐'!$ − 𝑉'!$,' 𝑦' 𝑦$ 𝐕 : , 1 + ⋯ + 𝑦'!% 𝐕 : , 𝑜 − 2 = 𝒄 −𝑦' 𝐕 : , 𝑜 − 𝑦'!$ 𝐕 : , 𝑜 − 1

Or in general (backward-substitution for upper triangular systems): 𝑦&= 𝑐& − ∑'(&)!

%

𝑉&'𝑦' 𝑉&& , 𝑗 = 𝑜 − 1, 𝑜 − 2, … , 1 𝑦% = 𝑐%/𝑉%%

Forward-substitution for lower-triangular systems: 𝑦&= 𝑐& − ∑'(!

&*! 𝑀&'𝑦'

𝑀&& , 𝑗 = 2,3, … , 𝑜 𝑦! = 𝑐!/𝑀!!

𝑀!! 𝑀"! 𝑀"" … … ⋮ ⋮ 𝑀#! 𝑀#" ⋱ ⋮ … 𝑀## 𝑦! 𝑦" ⋮ 𝑦# = 𝑐 𝑐" ⋮ 𝑐#

Triangular Matrices

Cost of solving triangular systems

𝑦&= 𝑐& − ∑'(&)!

%

𝑉&'𝑦' 𝑉&& , 𝑗 = 𝑜 − 1, 𝑜 − 2, … , 1 𝑦% = 𝑐%/𝑉%%

𝑜 divisions 𝑜 𝑜 − 1 /2 subtractions/additions 𝑜 𝑜 − 1 /2 multiplications

Computational complexity is 𝑃(𝑜")

𝑜 divisions 𝑜 𝑜 − 1 /2 subtractions/additions 𝑜 𝑜 − 1 /2 multiplications

Computational complexity is 𝑃(𝑜")

𝑦&= 𝑐& − ∑'(!

&*! 𝑀&'𝑦'

𝑀&& , 𝑗 = 2,3, … , 𝑜 𝑦! = 𝑐!/𝑀!!

Linear System of Equations

How do we solve 𝑩 𝒚 = 𝒄 when 𝑩 is a non-triangular matrix? We can perform LU factorization: given a 𝑜×𝑜 matrix 𝑩,

• btain lower triangular matrix 𝑴 and upper triangular matrix

𝑽 such that where we set the diagonal entries of 𝑴 to be equal to 1.

𝑩 = 𝑴𝑽

1 𝑀"! 1 … … ⋮ ⋮ 𝑀#! 𝑀#" ⋱ ⋮ … 1 𝑉!! 𝑉!" 𝑉"" … 𝑉!# … 𝑉"# ⋮ ⋮ ⋱ ⋮ … 𝑉## = 𝐵!! 𝐵!" 𝐵"! 𝐵"" … 𝐵!# … 𝐵"# ⋮ ⋮ 𝐵#! 𝐵#" ⋱ ⋮ … 𝐵##

LU Factorization

1 𝑀"! 1 … … ⋮ ⋮ 𝑀#! 𝑀#" ⋱ ⋮ … 1 𝑉!! 𝑉!" 𝑉"" … 𝑉!# … 𝑉"# ⋮ ⋮ ⋱ ⋮ … 𝑉## = 𝐵!! 𝐵!" 𝐵"! 𝐵"" … 𝐵!# … 𝐵"# ⋮ ⋮ 𝐵#! 𝐵#" ⋱ ⋮ … 𝐵##

𝑴𝑽 𝒚 = 𝒄 𝑴 𝒛 = 𝒄

Forward-substitution with complexity 𝑃(𝑜")

𝑽 𝒚 = 𝒛

Solve for 𝒛 Backward-substitution with complexity 𝑃(𝑜") Solve for 𝒚 Assuming the LU factorization is know, we can solve the general system By solving two triangular systems: But what is the cost of the LU factorization? Is it beneficial?

2x2 LU Factorization (simple example)

𝐵!! 𝐵!" 𝐵"! 𝐵"" = 1 𝑀"! 1 𝑉!! 𝑉!" 𝑉""

𝐵!! 𝐵!" 𝐵"! 𝐵"" = 𝑉!! 𝑉!" 𝑀"!𝑉!! 𝑀"!𝑉!" + 𝑉""

2) 𝑀=> = 𝐵=>/𝑉>> 3) 𝑉== = 𝐵== − 𝑀=>𝑉>= Seems quite simple! Can we generalize this for a 𝑜×𝑜 matrix 𝑩? 𝑉>> = 𝐵=>/𝑉>>

LU Factorization

𝐵!! 𝐵!" 𝐵"! 𝐵"" … 𝐵!# … 𝐵"# ⋮ ⋮ 𝐵#! 𝐵#" ⋱ ⋮ … 𝐵## = 𝑏!! 𝒃!" 𝒃"! 𝑩"" = 1 𝟏 𝒎"! 𝑴"" 𝑣!! 𝒗!" 𝟏 𝑽""

𝑏$$ 𝒃$% 𝒃%$ 𝑩%% 𝑏$$: scalar 𝒃$%: row vector (1×(𝑜 − 1)) 𝒃%$: column vector (𝑜 − 1)×1 𝑩%%: matrix (𝑜 − 1)×(𝑜 − 1)

𝑏!! 𝒃!" 𝒃"! 𝑩"" = 𝑣!! 𝒗!" 𝑣!! 𝒎"! 𝒎"!𝒗!" + 𝑴""𝑽""

1) First row of 𝑽 is the first row of 𝑩 𝒎=> = >

+!! 𝒃=>

3) 𝑴""𝑽"" = 𝑩"" − 𝒎"!𝒗!" Need another factorization!

Known!

2) First column of 𝑴 is the first column of 𝑩/ 𝑣>>

𝑩"" = 𝒎"!𝒗!" + 𝑴""𝑽""

LU Factorization

𝐵!! 𝐵!" 𝐵"! 𝐵"" … 𝐵!# … 𝐵"# ⋮ ⋮ 𝐵#! 𝐵#" ⋱ ⋮ … 𝐵## = 𝑏!! 𝒃!" 𝒃"! 𝑩"" = 1 𝟏 𝒎"! 𝑴"" 𝑣!! 𝒗!" 𝟏 𝑽""

𝑏$$ 𝒃$% 𝒃%$ 𝑩%% 𝑏$$: scalar 𝒃$%: row vector (1×(𝑜 − 1)) 𝒃%$: column vector (𝑜 − 1)×1 𝑩%%: matrix (𝑜 − 1)×(𝑜 − 1)

𝑏!! 𝒃!" 𝒃"! 𝑩"" = 𝑣!! 𝒗!" 𝑣!! 𝒎"! 𝒎"!𝒗!" + 𝑴""𝑽""

1) First row of 𝑽 is the first row of 𝑩 2) 𝒎=> = >

+!! 𝒃=>

3) 𝑵 = 𝑴""𝑽"" = 𝑩"" − 𝒎"!𝒗!" Need another factorization!

Known!

First column of 𝑴 is the first column of 𝑩/ 𝑣>>

Example

𝑵 = 2 8 1 2 4 1 3 3 1 2 1 3 6 2 4 2

1) First row of 𝑽 is the first row of 𝑩 2) First column of 𝑴 is the first column of 𝑩/ 𝑣>>

3) 𝑴""𝑽"" = 𝑩"" − 𝒎"!𝒗!"

Example

𝑵 = 2 8 1 2 4 1 3 3 1 2 1 3 6 2 4 2 𝑽 = 2 8 4 1 𝑴 = 1 0.5 0.5 0.5 𝑵 = 2 8 1 −2 4 1 1 2.5 1 −2 1 −1 4 1.5 2 1.5 𝑽 = 2 8 −2 4 1 1 2.5 𝑴 = 1 0.5 1 0.5 1 0.5 0.5 𝑵 = 2 8 1 −2 4 1 1 2.5 1 −2 1 −1 3 −1 1.5 0.25 𝑽 = 2 8 −2 4 1 1 2.5 3 −1 𝑴 = 1 0.5 1 0.5 1 0.5 0.5 1 0.5 𝑽 = 2 8 −2 4 1 1 2.5 3 −1 0.75 𝑴 = 1 0.5 1 0.5 1 0.5 0.5 1 0.5 1 𝑵 = 2 8 1 −2 4 1 1 2.5 1 −2 1 −1 3 −1 1.5 0.75

Algorithm: LU Factorization of matrix A

Cost of LU factorization

$

"#$ %

𝑗 = 1 2 𝑛 𝑛 + 1 $

"#$ %

𝑗& = 1 6 𝑛 𝑛 + 1 2𝑛 + 1

Side note:

Cost of LU factorization

Number of divisions: 𝑜 − 1 + 𝑜 − 2 + ⋯ + 1 = 𝑜 𝑜 − 1 /2 Number of multiplications 𝑜 − 1 % + 𝑜 − 2 % + … + 1 % =

'! & − '" % + ' "

Number of subtractions: 𝑜 − 1 % + 𝑜 − 2 % + … + 1 % =

'! & − '" % + ' "

Computational complexity is 𝑃(𝑜,)

$

"#$ %

𝑗 = 1 2 𝑛 𝑛 + 1 $

"#$ %

𝑗& = 1 6 𝑛 𝑛 + 1 2𝑛 + 1

Side note: Demo “Complexity of Mat-Mat multiplication and LU”

Solving linear systems

In general, we can solve a linear system of equations following the steps: 1) Factorize the matrix 𝑩 : 𝑩 = 𝑴𝑽 (complexity 𝑃(𝑜C))

2) Solve 𝑴 𝒛 = 𝒄 (complexity 𝑃(𝑜=))

3) Solve 𝑽 𝒚 = 𝒛 (complexity 𝑃(𝑜=))

But why should we decouple the factorization from the actual solve? (Remember from Linear Algebra, Gaussian Elimination does not decouple these two steps…)

Example: Optimization of automotive control arm

Find the distribution of material inside the design space (𝒆) that maximizes the stiffness, i.e., min 𝑽,𝑮 where 𝑳 𝒆 𝑽 = 𝑮 (𝑽: displacement vector, 𝑮: load vector, 𝑳: stiffness matrix) Solve the linear system of equations 𝑳 𝑽 = 𝑮 for the load vector 𝑮. What if we have many different loading conditions (pothole, hitting a curb, breaking, etc)?

Iclicker question

Let’s assume that when solving the system of equations 𝑳 𝑽 = 𝑮, we observe the following:

• When the stiffness matrix has dimensions (100,100), computing the LU factorization

takes about 1 second and each solve (forward + backward substitution) takes about 0.01 seconds. Estimate the total time it will take to find the displacement response corresponding to 10 different load vectors 𝑮 when the stiffness matrix has dimensions (1000,1000)? 𝐵) ~10 𝑡𝑓𝑑𝑝𝑜𝑒𝑡 𝐶) ~10" 𝑡𝑓𝑑𝑝𝑜𝑒𝑡 𝐷) ~10# 𝑡𝑓𝑑𝑝𝑜𝑒𝑡 𝐸) ~10$ 𝑡𝑓𝑑𝑝𝑜𝑒𝑡 𝐹) ~10- 𝑡𝑓𝑑𝑝𝑜𝑒𝑡

What can go wrong with the previous algorithm?

If division by zero occurs, LU factorization fails. What can we do to get something like an LU factorization?

What can go wrong with the previous algorithm?

𝑵 = 2 8 1 𝟓 4 1 3 3 1 2 1 3 6 2 4 2 𝑽 = 2 8 4 1 𝑴 = 1 0.5 0.5 0.5 𝑵 − 𝒎"!𝒗!" = 2 8 1 𝟏 4 1 1 2.5 1 −2 1 −1 4 1.5 2 1.5 𝒎"!𝒗!" = 4 2 0.5 4 2 0.5 4 2 0.5

The next update for the lower triangular matrix will result in a division by zero! LU factorization fails. What can we do to get something like an LU factorization?

Demo “Little c”

Pivoting

Approach:

• 1. Swap rows if there is a zero entry in the diagonal
• 2. Even better idea: Find the largest entry (by absolute value) and

swap it to the top row. The entry we divide by is called the pivot. Swapping rows to get a bigger pivot is called (partial) pivoting.

𝑏!! 𝒃!" 𝒃"! 𝑩"" = 𝑣!! 𝒗!" 𝑣!! 𝒎"! 𝒎"!𝒗!" + 𝑴""𝑽"" Find the largest entry (in magnitude)

LU Factorization with Partial Pivoting

• LU factorization with partial pivoting can be completed for any matrix A

Suppose you are at stage k and there is no non-zero entry on or below the diagonal in column k. At this point, there is nothing else you can do, so the algorithm leaves a zero in the diagonal entry of U. Note that the matrix U is singular, and so is the matrix A. Subsequent backward substitutions using U will fail, but the LU factorization itself is still completed.

LU Factorization with Partial Pivoting

where 𝑸 is a permutation matrix Then solve two triangular systems:

𝑩 = 𝑸𝑴𝑽 𝑴 𝒛 = 𝑸#𝒄 𝑽 𝒚 = 𝒛

(Solve for 𝒛) (Solve for 𝒚)

𝑩 𝒚 = 𝒄 → 𝑸𝑴𝑽 𝒚 = 𝒄 → 𝑴𝑽 𝒚 = 𝑸#𝒄

Example

𝑩 = 𝑵 = 2 8 1 2 4 1 3 3 1 2 1 3 3 2 4 2 𝑽 = 2 8 4 1 𝑴 = 1 0.5 0.5 0.5 𝑵 = 2 8 1 −2 4 1 1 2.5 1 −2 1 −1 1 1.5 2 1.5 𝑽 = 2 8 −2 4 1 1 2.5 𝑴 = 1 0.5 1 0.5 1 0.5 0.5 𝑵 = 2 8 1 −2 4 1 1 2.5 1 −2 1 −1 −1 1.5 0.25 𝑵 = 2 8 1 −2 4 1 1 2.5 1 −1 1 −2 1.5 0.25 −1 𝑽 = 2 8 −2 4 1 1 2.5 1.5 0.25 −1 𝑽 = 2 8 −2 4 1 1 2.5 𝑴 = 1 0.5 1 0.5 0.5 0.5 1.0 𝑴 = 1 0.5 1 0.5 0.5 0.5 1.0 1.0 1.0

Demo “Pivoting example”

𝑩 = 2 1 4 3 1 3 1 8 7 6 7 9 5 9 8 𝑽 = 8 7 9 5 𝑴 = 1 0.5 0.25 0.75 Y 𝑩 = 𝑸𝑩 = 1 1 1 1 2 1 4 3 1 3 1 8 7 6 7 9 5 9 8 = 8 7 4 3 9 5 3 1 2 1 6 7 1 9 8

D 𝑩 − 𝒎%$𝒗$% = 8 7 4 −0.5 9 5 −1.5 −1.5 2 −0.75 6 1.75 −1.25 −1.25 2.25 4.25 𝒎%$𝒗$% = 3.5 4.5 2.5 1.75 2.25 1.25 5.25 6.75 3.75

Demo “Pivoting example”

𝑽 = 8 7 1.75 9 5 2.25 4.25 𝑴 = 1 0.75 1 0.25 −0.428 0.5 −0.285 7 𝑩 = 𝑸7 𝑩 = 1 1 1 1 8 7 6 1.75 9 5 2.25 4.25 2 −0.75 4 −0.5 −1.25 −1.25 −1.5 −1.5 = 8 7 6 1.75 9 5 2.25 4.25 2 −0.75 4 −0.5 −1.25 −1.25 −1.5 −1.5 7 𝑩 = 7 𝑩 − 𝒎&$𝒗$& = 8 7 4 −0.5 9 5 −1.5 −1.5 2 −0.75 6 1.75 −1.25 −1.25 2.25 4.25 𝒎&$𝒗$& = −0.963 −1.819 −0.6412 −1.2112 7 𝑩 = 7 𝑩 − 𝒎&$𝒗$& = 8 7 6 1.75 9 5 2.25 4.25 2 −0.75 4 −0.5 −0.287 0.569 −0.8587 −0.2887

Demo “Pivoting example”

𝑽 = 8 7 1.75 9 5 2.25 4.25 −0.86 −0.29 𝑴 = 1 0.75 1 0.5 −0.285 0.25 −0.428 1 0.334 7 𝑩 = 𝑸7 𝑩 = 1 1 1 1 8 7 6 1.75 9 5 2.25 4.25 2 −0.75 4 −0.5 −0.287 0.569 −0.8587 −0.2887 = 8 7 6 1.75 9 5 2.25 4.25 4 −0.5 2 −0.75 −0.8587 −0.2887 −0.287 0.569 7 𝑩 = 7 𝑩 − 𝒎&$𝒗$& = 8 7 6 1.75 9 5 2.25 4.25 2 −0.75 4 −0.5 −0.287 0.569 −0.8587 −0.2887 𝑽 = 8 7 1.75 9 5 2.25 4.25 −0.86 −0.29 0.67 𝑴 = 1 0.75 1 0.5 −0.285 0.25 −0.428 1 0.334 1

𝑸 = 1 1 1 1