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” What happens if we add some noise to 𝒄? 𝒚 𝒄

”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

Example: Backward-substitution for upper

triangular systems

2 8 4 4 2 4 3 6 2 2 𝑦! 𝑦" 𝑦$ 𝑦% = 2 4 4 1

𝑦# = 1 2 𝑦& = 4 − 2 1 2 6 = 1 2 𝑦% = 4 − 4 1 2 − 3 1 2 4 = 1/2 4 = 1 8 𝑦$ = 2 − 8 1 8 − 4 1 2 − 2 1 2 2 = −2 2 = −1

LU Factorization

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 𝑉!! 𝑉!" 𝑉"" … 𝑉!# … 𝑉"# ⋮ ⋮ ⋱ ⋮ … 𝑉## = 𝐵!! 𝐵!" 𝐵"! 𝐵"" … 𝐵!# … 𝐵"# ⋮ ⋮ 𝐵#! 𝐵#" ⋱ ⋮ … 𝐵## Assuming the LU factorization is know, we can solve the general system

LU Factorization (with pivoting)

𝑸𝑴𝑽 𝒚 = 𝒄 𝑴 𝒛 = 𝑸𝑼𝒄

Forward-substitution

𝑽 𝒚 = 𝒛

(Solve for 𝒛) Backward-substitution (Solve for 𝒚)

𝑩 = 𝑸𝑴𝑽

Factorize:

𝒛

Example

𝑽 = 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

Assume the 𝑩 = 𝑴𝑽 factorization is known, yielding: Determine the solution 𝒚 that satisfies 𝑩𝒚 = 𝒄, when 𝒄 = 2 2 1 4 First, solve the lower-triangular system 𝑴 𝒛 = 𝒄 for the variable 𝒛 𝑴 𝑽𝒚 = 𝒄

𝒛

Then, solve the upper-triangular system 𝑽 𝒚 = 𝒛 for the variable 𝒚

1 0.5 1 0.5 1 0.5 0.5 1 0.5 1 𝒛 = 2 2 1 4 2 8 −2 4 1 1 2.5 3 −1 0.75 𝒚 = 2 1 −1 3

Methods to solve linear system of equations

𝑩 𝒚 = 𝒄

• LU
• Cholesky
• Sparse

LU Factorization - Algorithm

2x2 LU Factorization (simple example)

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

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

LU Factorization

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

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

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

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 𝑽 = 2 8 4 1 𝑴 = 1 0.5 0.5 0.5

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

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

𝑴%%𝑽%% = 𝑩%% − 𝒎%$𝒗$% = 2 3 3 2 3 6 2 4 2 − 4 2 0.5 4 4 2 0.5 2 0.5

𝑵 = 2 8 1 −2 4 1 1 2.5 1 −2 1 −1 4 1.5 2 1.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

𝑴%%𝑽%% = 𝑩%% − 𝒎%$𝒗$% = 4 1.5 2 1.5 − 1 2.5 0.5 1.25

𝑵 = 2 8 1 −2 4 1 1 2.5 1 −2 1 −1 3 −1 1.5 0.25

𝑵 = 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

𝑴%%𝑽%% = 𝑩%% − 𝒎%$𝒗$% = 0.25 − −0.5 = 0.75

𝑽 = 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

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 𝑩/ 𝑣==

Cost of solving linear system of equations

Cost of solving triangular systems

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

*

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

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, … , 𝑜 𝑦! = 𝑐!/𝑀!!

Cost of LU factorization

&

"#$ %

𝑗 = 1 2 𝑛 𝑛 + 1 &

"#$ %

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

Side note:

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

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

• When the 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 response 𝑽 corresponding to 10 different vectors 𝑮 when the matrix 𝑳 has dimensions (1000,1000)? 𝐵) ~10 𝑡𝑓𝑑𝑝𝑜𝑒𝑡 𝐶) ~10" 𝑡𝑓𝑑𝑝𝑜𝑒𝑡 𝐷) ~10# 𝑡𝑓𝑑𝑝𝑜𝑒𝑡 𝐸) ~10$ 𝑡𝑓𝑑𝑝𝑜𝑒𝑡 𝐹) ~10, 𝑡𝑓𝑑𝑝𝑜𝑒𝑡

LU Factorization with pivoting

What can go wrong with the previous algorithm for LU factorization?

𝑵 = 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?

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)