Numerical Linear Algebra EECS 442 David Fouhey Fall 2019, - PowerPoint PPT Presentation

Matrices Horizontally concatenate n, m-dim column vectors and you get a mxn matrix A (here 2x3) 𝑤 1 1 𝑤 2 1 𝑤 3 1 𝑩 = 𝒘 1 , ⋯ , 𝒘 𝑜 = 𝑤 1 2 𝑤 2 2 𝑤 3 2 (scalar) (vector) (matrix) a undecorated a bold or arrow A lowercase lowercase uppercase bold

Matrices 𝑈 𝑏 Transpose: flip (3x1) T = 1x3 𝑐 = 𝑏 𝑐 𝑑 rows / columns 𝑑 Vertically concatenate m, n-dim row vectors and you get a mxn matrix A (here 2x3) 𝑈 𝒗 1 𝑣 1 1 𝑣 1 2 𝑣 1 3 𝐵 = = ⋮ 𝑣 2 1 𝑣 2 2 𝑣 2 3 𝑈 𝒗 𝑜

Matrix-Vector Product 𝒛 2𝑦1 = 𝑩 2𝑦3 𝒚 3𝑦1 𝑦 1 𝑧 1 𝑦 2 𝒘 𝟐 𝒘 𝟑 𝒘 𝟒 𝑧 2 = 𝑦 3 𝒛 = 𝑦 1 𝒘 𝟐 + 𝑦 2 𝒘 𝟑 + 𝑦 3 𝒘 𝟒 Linear combination of columns of A

Matrix-Vector Product 𝒛 2𝑦1 = 𝑩 2𝑦3 𝒚 3𝑦1 3 𝑧 1 𝒗 𝟐 𝑧 2 = 𝒚 3 𝒗 𝟑 𝑼 𝒚 𝑼 𝒚 𝑧 1 = 𝒗 𝟐 𝑧 2 = 𝒗 𝟑 Dot product between rows of A and x

Matrix Multiplication Generally: A mn and B np yield product ( AB ) mp | | 𝑼 − 𝒃 𝟐 − 𝒄 𝟐 ⋯ 𝒄 𝒒 𝑩𝑪 = ⋮ 𝑼 − 𝒃 𝒏 − | | Yes – in A , I’m referring to the rows, and in B , I’m referring to the columns

Matrix Multiplication Generally: A mn and B np yield product ( AB ) mp | | 𝒄 𝟐 ⋯ 𝒄 𝒒 | | 𝑼 𝒄 𝟐 𝑼 𝒄 𝒒 𝑼 𝒃 𝟐 ⋯ 𝒃 𝟐 − 𝒃 𝟐 − 𝑩𝑪 = ⋮ ⋮ ⋱ ⋮ 𝑼 𝒄 𝟐 𝑼 𝒄 𝒒 𝑼 − 𝒃 𝒏 − 𝒃 𝒏 ⋯ 𝒃 𝒏 𝑼 𝒄 𝒌 𝑩𝑪 𝑗𝑘 = 𝒃 𝒋

Matrix Multiplication • Dimensions must match • Dimensions must match • Dimensions must match • (Yes, it’s associative): ABx = (A)(Bx) = (AB)x • ( No it’s not commutative ): ABx ≠ (BA)x ≠ ( BxA)

Operations They Don’t Teach You Probably Saw Matrix Addition 𝑒 + 𝑓 𝑔 𝑏 + 𝑓 𝑐 + 𝑔 𝑏 𝑐 ℎ = 𝑕 𝑑 + 𝑕 𝑒 + ℎ 𝑑 What is this? FYI: e is a scalar 𝑏 𝑐 𝑏 + 𝑓 𝑐 + 𝑓 𝑒 + 𝑓 = 𝑑 𝑑 + 𝑓 𝑒 + 𝑓

Broadcasting If you want to be pedantic and proper, you expand e by multiplying a matrix of 1s (denoted 1 ) 𝑏 𝑐 = 𝑏 𝑐 𝑒 + 𝑓 𝑒 + 𝟐 2𝑦2 𝑓 𝑑 𝑑 𝑒 + 𝑓 𝑓 = 𝑏 𝑐 𝑓 𝑓 𝑑 Many smart matrix libraries do this automatically. This is the source of many bugs.

Broadcasting Example Given: a nx2 matrix P and a 2D column vector v , Want: nx2 difference matrix D 𝑦 1 𝑧 1 𝑦 1 − 𝑏 𝑧 1 − 𝑐 𝒘 = 𝑏 ⋮ ⋮ 𝑸 = ⋮ ⋮ 𝑬 = 𝑐 𝑦 𝑜 𝑧 𝑜 𝑦 𝑜 − 𝑏 𝑧 𝑜 − 𝑐 𝑦 1 𝑧 1 Blue stuff is 𝑏 𝑐 𝑸 − 𝒘 𝑈 = ⋮ ⋮ assumed / − ⋮ broadcast 𝑦 𝑜 𝑧 𝑜 𝑏 𝑐

Two Uses for Matrices 1. Storing things in a rectangular array (images, maps) • Typical operations : element-wise operations, convolution (which we’ll cover next) • Atypical operations : almost anything you learned in a math linear algebra class 2. A linear operator that maps vectors to another space ( Ax ) • Typical/Atypical: reverse of above

Images as Matrices Suppose someone hands you this matrix. What’s wrong with it?

Contrast – Gamma curve Typical way to change the contrast is to apply a nonlinear correction pixelvalue 𝛿 The quantity 𝛿 controls how much contrast gets added

Contrast – Gamma curve Now the darkest 90% regions (10 th pctile) are much darker 50% than the moderately new dark regions (50 th 10% 90% pctile). new new 10% 50%

Implementation Python+Numpy (right way): imNew = im**0.25 Python+Numpy (slow way – why? ): imNew = np.zeros(im.shape) for y in range(im.shape[0]): for x in range(im.shape[1]): imNew[y,x] = im[y,x]**expFactor

Results Phew! Much Better.

Element-wise Operations Element-wise power – beware notation 𝑞 𝑩 𝑞 𝑗𝑘 = 𝐵 𝑗𝑘 “Hadamard Product” / Element -wise multiplication 𝑩 ⊙ 𝑪 𝑗𝑘 = 𝑩 𝑗𝑘 ∗ 𝑪 𝑗𝑘 Element-wise division 𝑩/𝑪 𝑗𝑘 = 𝐵 𝑗𝑘 𝐶 𝑗𝑘

Sums Across Axes 𝑦 1 𝑧 1 Suppose have ⋮ ⋮ 𝑩 = Nx2 matrix A 𝑦 𝑜 𝑧 𝑜 𝑦 1 + 𝑧 1 ND col. vec. ⋮ Σ(𝑩, 1) = 𝑦 𝑜 + 𝑧 𝑜 𝑜 𝑜 2D row vec Σ(𝑩, 0) = ෍ 𝑦 𝑗 , ෍ 𝑧 𝑗 𝑗=1 𝑗=1 Note – libraries distinguish between N-D column vector and Nx1 matrix.

Vectorizing Example • Suppose I represent each image as a 128- dimensional vector • I want to compute all the pairwise distances between { x 1 , …, x N } and { y 1 , …, y M } so I can find, for every x i the nearest y j • Identity: 𝒚 − 𝒛 2 = 𝒚 2 + 𝒛 2 − 2𝒚 𝑈 𝒛 𝒚 2 + 𝒛 2 − 2𝒚 𝑈 𝒛 1/2 • Or: 𝒚 − 𝒛 =

Vectorizing Example − 𝒚 1 − − 𝒛 1 − | | 𝒁 𝑼 = ⋮ ⋮ 𝒛 1 ⋯ 𝒛 𝑁 𝒀 = 𝒁 = − 𝒚 𝑂 − − 𝒛 𝑁 − | | 𝒚 𝟐 𝟑 Compute a Nx1 𝚻 𝒀 𝟑 , 𝟐 = ⋮ vector of norms 𝒚 𝑶 𝟑 (can also do Mx1) Compute a NxM 𝑼 𝒛 𝒌 𝒀𝒁 𝑼 𝑗𝑘 = 𝒚 𝒋 matrix of dot products

Vectorizing Example 𝑼 − 2𝒀𝒁 𝑼 1/2 𝐄 = Σ 𝒀 𝟑 , 1 + Σ 𝒁 𝟑 , 1 𝒚 𝟐 𝟑 𝒛 1 𝟑 𝟑 + ⋯ 𝒛 𝑁 ⋮ 𝒚 𝑶 𝟑 𝒚 𝟐 2 + 𝒛 𝟐 2 𝒚 𝟐 2 + 𝒛 𝑵 2 ⋯ Why? ⋮ ⋱ ⋮ 𝒚 𝑶 2 + 𝒛 𝟐 2 𝒚 𝑶 2 + 𝒛 𝑵 2 ⋯ 2 2 + 𝑧 𝑘 Σ 𝑌 2 , 1 + Σ 𝑍 2 , 1 𝑈 𝑗𝑘 = 𝑦 𝑗

Vectorizing Example 𝑼 − 2𝒀𝒁 𝑼 1/2 𝐄 = Σ 𝒀 𝟑 , 1 + Σ 𝒁 𝟑 , 1 2 + 2𝒚 𝑼 𝒛 𝒚 𝒋 2 + 𝒛 𝒌 𝐄 𝑗𝑘 = Numpy code: XNorm = np.sum(X**2,axis=1,keepdims=True) YNorm = np.sum(Y**2,axis=1,keepdims=True) D = (XNorm+YNorm.T-2*np.dot(X,Y.T))**0.5 *May have to make sure this is at least 0 (sometimes roundoff issues happen)

Does it Make a Difference? Computing pairwise distances between 300 and 400 128-dimensional vectors 1. for x in X, for y in Y, using native python: 9s 2. for x in X, for y in Y, using numpy to compute distance: 0.8s 3. vectorized: 0.0045s (~2000x faster than 1, 175x faster than 2) Expressing things in primitives that are optimized is usually faster

Linear Independence A set of vectors is linearly independent if you can’t write one as a linear combination of the others. 0 0 5 Suppose: 𝒃 = 𝒄 = 𝒅 = 0 6 0 2 0 0 0 0 = 1 2 𝒃 − 1 𝒚 = = 2𝒃 0 𝒛 = 3 𝒄 −2 4 1 • Is the set {a,b,c} linearly independent? • Is the set {a,b,x} linearly independent? • Max # of independent 3D vectors?

Span Span: all linear combinations of a set of vectors Span({ }) = Span({[0,1]}) = ? All vertical lines through origin = 𝜇 0,1 : 𝜇 ∈ 𝑆 Is blue in {red }’s span?

Span Span: all linear combinations of a set of vectors Span({ , }) = ?

Span Span: all linear combinations of a set of vectors Span({ , }) = ?

Matrix-Vector Product | | Right-multiplying A by x 𝒅 𝟐 ⋯ 𝒅 𝒐 mixes columns of A 𝑩𝒚 = 𝒚 according to entries of x | | • The output space of f( x ) = Ax is constrained to be the span of the columns of A . • Can’t output things you can’t construct out of your columns

An Intuition 𝑦 1 | | | 𝑦 2 𝒅 𝟐 𝒅 𝟑 𝒅 𝒐 𝒛 = 𝑩𝒚 = 𝑦 3 | | | y 1 y x Ax y 2 x 1 x 2 x 3 y 3 x – knobs on machine (e.g., fuel, brakes) y – state of the world (e.g., where you are) A – machine (e.g., your car)

Linear Independence Suppose the columns of 3x3 matrix A are not linearly independent (c 1 , α c 1 , c 2 for instance) 𝑦 1 | | | 𝑦 2 𝒅 𝟐 𝛽𝒅 𝟐 𝒅 𝟑 𝒛 = 𝑩𝒚 = 𝑦 3 | | | 𝒛 = 𝑦 1 𝒅 𝟐 + 𝛽𝑦 2 𝒅 𝟐 + 𝑦 3 𝒅 𝟑 𝒛 = 𝑦 1 + 𝛽𝑦 2 𝒅 𝟐 + 𝑦 3 𝒅 𝟑

Linear Independence Intuition Knobs of x are redundant. Even if y has 3 outputs, you can only control it in two directions 𝒛 = 𝑦 1 + 𝛽𝑦 2 𝒅 𝟐 + 𝑦 3 𝒅 𝟑 y 1 y x Ax y 2 x 1 x 2 x 3 y 3

Linear Independence Recall: 𝑩𝒚 = 𝑦 1 + 𝛽𝑦 2 𝒅 𝟐 + 𝑦 3 𝒅 𝟑 𝑦 1 + 𝛾 𝑦 1 + 𝛾 + 𝛽𝑦 2 − 𝛽 𝛾 𝒛 = 𝑩 𝑦 2 − 𝛾/𝛽 = 𝛽 𝑑 1 + 𝑦 3 𝑑 2 𝑦 3 • Can write y an infinite number of ways by 𝛾 adding 𝛾 to x 1 and subtracting 𝛽 from x 2 • Or, given a vector y there’s not a unique vector x s.t. y = Ax • Not all y have a corresponding x s.t. y=Ax

Linear Independence 𝑩𝒚 = 𝑦 1 + 𝛽𝑦 2 𝒅 𝟐 + 𝑦 3 𝒅 𝟑 𝛾 𝛾 − 𝛽 𝛾 𝒛 = 𝑩 −𝛾/𝛽 = 𝛽 𝒅 𝟐 + 0𝒅 𝟑 0 • What else can we cancel out? • An infinite number of non-zero vectors x can map to a zero-vector y • Called the right null-space of A.

Rank • Rank of a nxn matrix A – number of linearly independent columns ( or rows ) of A / the dimension of the span of the columns • Matrices with full rank (n x n, rank n) behave nicely: can be inverted, span the full output space, are one-to-one. • Matrices with full rank are machines where every knob is useful and every output state can be made by the machine

Inverses • Given 𝒛 = 𝑩𝒚 , y is a linear combination of columns of A proportional to x . If A is full-rank, we should be able to invert this mapping. • Given some y (output) and A , what x (inputs) produced it? • x = A -1 y • Note: if you don’t need to compute it, never ever compute it. Solving for x is much faster and stable than obtaining A -1 .

Symmetric Matrices • Symmetric: 𝑩 𝑼 = 𝑩 or 𝑏 11 𝑏 12 𝑏 13 𝑩 𝑗𝑘 = 𝑩 𝑘𝑗 𝑏 21 𝑏 22 𝑏 23 • Have lots of special 𝑏 31 𝑏 32 𝑏 33 properties Any matrix of the form 𝑩 = 𝒀 𝑼 𝒀 is symmetric. 𝑼 𝑩 𝑼 = 𝒀 𝑼 𝒀 Quick check: 𝑩 𝑼 = 𝒀 𝑼 𝒀 𝑼 𝑼 𝑩 𝑼 = 𝒀 𝑼 𝒀

Special Matrices – Rotations 𝑠 𝑠 𝑠 11 12 13 𝑠 𝑠 𝑠 21 22 23 𝑠 𝑠 𝑠 31 32 33 • Rotation matrices 𝑺 rotate vectors and do not change vector L2 norms ( 𝑺𝒚 2 = 𝒚 2 ) • Every row/column is unit norm • Every row is linearly independent • Transpose is inverse 𝑺𝑺 𝑼 = 𝑺 𝑼 𝑺 = 𝑱 • Determinant is 1 (otherwise it’s also a coordinate flip/reflection), eigenvalues are 1

Eigensystems • An eigenvector 𝒘 𝒋 and eigenvalue 𝜇 𝑗 of a matrix 𝑩 satisfy 𝑩𝒘 𝒋 = 𝜇 𝑗 𝒘 𝒋 ( 𝑩𝒘 𝒋 is scaled by 𝜇 𝑗 ) • Vectors and values are always paired and typically you assume 𝒘 𝒋 2 = 1 • Biggest eigenvalue of A gives bounds on how much 𝑔 𝒚 = 𝑩𝒚 stretches a vector x . • Hints of what people really mean: • “Largest eigenvector” = vector w/ largest value • Spectral just means there’s eigenvectors

Suppose I have points in a grid

Now I apply f( x ) = Ax to these points Pointy-end: Ax . Non-Pointy-End: x

𝑩 = 1.1 0 0 1.1 Red box – unit square, Blue box – after f( x ) = Ax . What are the yellow lines and why?

𝑩 = 0.8 0 0 1.25 Now I apply f( x ) = Ax to these points Pointy-end: Ax . Non-Pointy-End: x

𝑩 = 0.8 0 0 1.25 Red box – unit square, Blue box – after f( x ) = Ax . What are the yellow lines and why?

𝑩 = cos(𝑢) −sin(𝑢) sin(𝑢) cos(𝑢) Red box – unit square, Blue box – after f( x ) = Ax . Can we draw any yellow lines?

Eigenvectors of Symmetric Matrices • Always n mutually orthogonal eigenvectors with n (not necessarily) distinct eigenvalues • For symmetric 𝑩 , the eigenvector with the 𝒚 𝑼 𝑩𝒚 largest eigenvalue maximizes 𝒚 𝑼 𝒚 (smallest/min) • So for unit vectors (where 𝒚 𝑼 𝒚 = 1 ), that eigenvector maximizes 𝒚 𝑼 𝑩𝒚 • A surprisingly large number of optimization problems rely on (max/min)imizing this

The Singular Value Decomposition Can always write a mxn matrix A as: 𝑩 = 𝑽𝚻𝑾 𝑼 0 σ 1 A = U ∑ σ 2 σ 3 Scale 0 Rotation Eigenvectors Sqrt of of AA T Eigenvalues of A T A

The Singular Value Decomposition Can always write a mxn matrix A as: 𝑩 = 𝑽𝚻𝑾 𝑼 V T A = U ∑ Rotation Scale Rotation Eigenvectors Sqrt of Eigenvectors of AA T Eigenvalues of A T A of A T A

Singular Value Decomposition • Every matrix is a rotation, scaling, and rotation • Number of non-zero singular values = rank / number of linearly independent vectors • “Closest” matrix to A with a lower rank 0 σ 1 V T σ 2 = A U σ 3 0

Singular Value Decomposition • Every matrix is a rotation, scaling, and rotation • Number of non-zero singular values = rank / number of linearly independent vectors • “Closest” matrix to A with a lower rank 0 σ 1 V T σ 2 = Â U 0 0

Singular Value Decomposition • Every matrix is a rotation, scaling, and rotation • Number of non-zero singular values = rank / number of linearly independent vectors • “Closest” matrix to A with a lower rank • Secretly behind basically many things you do with matrices

Solving Least-Squares Start with two points (x i ,y i ) (x 2 ,y 2 ) 𝒛 = 𝑩𝒘 𝑧 1 𝑧 2 = 𝑦 1 1 𝑛 𝑦 2 1 𝑐 (x 1 ,y 1 ) 𝑧 1 𝑧 2 = 𝑛𝑦 1 + 𝑐 𝑛𝑦 2 + 𝑐 We know how to solve this – invert A and find v (i.e., (m,b) that fits points)

Solving Least-Squares Start with two points (x i ,y i ) (x 2 ,y 2 ) 𝒛 = 𝑩𝒘 𝑧 1 𝑧 2 = 𝑦 1 1 𝑛 𝑦 2 1 𝑐 (x 1 ,y 1 ) 2 𝑧 1 𝑧 2 − 𝑛𝑦 1 + 𝑐 𝒛 − 𝑩𝒘 2 = 𝑛𝑦 2 + 𝑐 2 + 𝑧 2 − 𝑛𝑦 2 + 𝑐 2 = 𝑧 1 − 𝑛𝑦 1 + 𝑐 The sum of squared differences between the actual value of y and what the model says y should be.

Solving Least-Squares Suppose there are n > 2 points 𝒛 = 𝑩𝒘 𝑧 1 𝑦 1 1 𝑛 ⋮ ⋮ ⋮ = 𝑐 𝑧 𝑂 𝑦 𝑂 1 Compute 𝑧 − 𝐵𝑦 2 again 𝑜 𝒛 − 𝑩𝒘 2 = ෍ 𝑧 𝑗 − (𝑛𝑦 𝑗 + 𝑐) 2 𝑗=1

Solving Least-Squares Given y , A , and v with y = Av overdetermined ( A tall / more equations than unknowns) We want to minimize 𝒛 − 𝑩𝒘 𝟑 , or find: arg min 𝒘 𝒛 − 𝑩𝒘 2 (The value of x that makes the expression smallest) Solution satisfies 𝑩 𝑼 𝑩 𝒘 ∗ = 𝑩 𝑼 𝒛 or −1 𝑩 𝑼 𝒛 𝒘 ∗ = 𝑩 𝑼 𝑩 (Don’t actually compute the inverse!)

When is Least-Squares Possible? Given y , A , and v . Want y = Av Want n outputs, have n knobs y = A v to fiddle with, every knob is useful if A is full rank. A: rows (outputs) > columns = (knobs). Thus can’t get precise v y A output you want (not enough knobs). So settle for “closest” knob setting.

When is Least-Squares Possible? Given y , A , and v . Want y = Av Want n outputs, have n knobs y = A v to fiddle with, every knob is useful if A is full rank. A: columns (knobs) > rows y = A (outputs). Thus, any output can v be expressed in infinite ways.

Homogeneous Least-Squares Given a set of unit vectors (aka directions) 𝒚 𝟐 , … , 𝒚 𝒐 and I want vector 𝒘 that is as orthogonal to all the 𝒚 𝒋 as possible (for some definition of orthogonal) Stack 𝒚 𝒋 into A , compute Av 𝒚 𝒐 … 𝒚 𝟑 0 if 𝑼 𝑼 𝒘 − 𝒚 𝟐 − 𝒚 𝟐 𝒚 𝟐 orthog 𝑩𝒘 = 𝒘 = ⋮ ⋮ 𝑼 𝑼 𝒘 − 𝒚 𝒐 − 𝒚 𝒐 𝒐 𝒘 𝟑 𝑩𝒘 𝟑 = ෍ 𝑼 𝒘 Compute 𝒚 𝒋 𝒋 Sum of how orthog. v is to each x

Homogeneous Least-Squares • A lot of times, given a matrix A we want to find the v that minimizes 𝑩𝒘 2 . • I.e., want 𝐰 ∗ = arg min 2 𝑩𝒘 2 𝒘 • What’s a trivial solution? • Set v = 0 → Av = 0 • Exclude this by forcing v to have unit norm

Homogeneous Least-Squares 2 Let’s look at 𝑩𝒘 2 2 = Rewrite as dot product 𝑩𝒘 2 2 = 𝐁𝐰 T (𝐁𝐰) Distribute transpose 𝑩𝒘 2 2 = 𝒘 𝑼 𝑩 𝑼 𝐁𝐰 = 𝐰 𝐔 𝐁 𝐔 𝐁 𝐰 𝑩𝒘 2 We want the vector minimizing this quadratic form Where have we seen this?