You can read an overview of this Numerical Linear Algebra course in this blog post. The course was originally taught in the University of San Francisco MS in Analytics graduate program. Course lecture videos are available on YouTube (note that the notebook numbers and video numbers do not line up, since some notebooks took longer than 1 video to cover).

You can ask questions about the course on our fast.ai forums.

3. Background Removal with Robust PCA¶

Our goal today:

Getting Started¶

Let's use the real video 003 dataset from BMC 2012 Background Models Challenge Dataset

Other sources of datasets:

• Human Activity Video Datasets
• Background Subtraction Website (a few links on this site are broken/outdated, but many work)

Import needed libraries:

Helper Methods¶

Load and view the data¶

An image from 1 moment in time is 60 pixels by 80 pixels (when scaled). We can unroll that picture into a single tall column. So instead of having a 2D picture that is $60 \times 80$, we have a $1 \times 4,800$ column

This isn't very human-readable, but it's handy because it lets us stack the images from different times on top of one another, to put a video all into 1 matrix. If we took the video image every tenth of a second for 113 seconds (so 11,300 different images, each from a different point in time), we'd have a $11300 \times 4800$ matrix, representing the video!

Since create_data_from_matrix is somewhat slow, we will save our matrix. In general, whenever you have slow pre-processing steps, it's a good idea to save the results for future use.

Note: High-res M is too big to plot, so only run the below with the low-res version

Questions: What are those wavy black lines? What are the horizontal lines?

SVD¶

Review from Lesson 2:

• What kind of matrices are returned by SVD?
• What is a way to speed up truncated SVD?

A first attempt with SVD¶

The below images were created with high-res data. Very slow to process:

Rank 1 Approximation¶

Let's zoom in on the people:

This is amazing for such a simple approach! We get somewhat better results through a more complicated algorithm below.

Principal Component Analysis (PCA)¶

When dealing with high-dimensional data sets, we often leverage on the fact that the data has low intrinsic dimensionality in order to alleviate the curse of dimensionality and scale (perhaps it lies in a low-dimensional subspace or lies on a low-dimensional manifold). Principal component analysis is handy for eliminating dimensions. Classical PCA seeks the best rank-$k$ estimate $L$ of $M$ (minimizing $\| M - L \|$ where $L$ has rank-$k$). Truncated SVD makes this calculation!

Traditional PCA can handle small noise, but is brittle with respect to grossly corrupted observations-- even one grossly corrupt observation can significantly mess up answer.

Robust PCA factors a matrix into the sum of two matrices, $M = L + S$, where $M$ is the original matrix, $L$ is low-rank, and $S$ is sparse. This is what we'll be using for the background removal problem! Low-rank means that the matrix has a lot of redundant information-- in this case, it's the background, which is the same in every scene (talk about redundant info!). Sparse means that the matrix has mostly zero entries-- in this case, see how the picture of the foreground (the people) is mostly empty. (In the case of corrupted data, $S$ is capturing the corrupted entries).

Applications of Robust PCA¶

• Video Surveillance

• Face Recognition photos from this excellent tutorial. The dataset here consists of images of the faces of several people taken from the same angle but with different illumniations.

  ([Source: Jean Kossaifi](https://jeankossaifi.github.io/tensorly/rpca.html)) ([Source: Jean Kossaifi](https://jeankossaifi.github.io/tensorly/rpca.html))

• Latent Semantic Indexing: $L$ captures common words used in all documents while $S$ captures the few key words that best distinguish each document from others

• Ranking and Collaborative Filtering: a small portion of the available rankings could be noisy and even tampered with (see Netflix RAD - Outlier Detection on Big Data on the official netflix blog)

L1 norm induces sparsity¶

The unit ball $\lVert x \rVert_1 = 1$ is a diamond in the L1 norm. It's extrema are the corners:

([Source](https://www.quora.com/Why-is-L1-regularization-supposed-to-lead-to-sparsity-than-L2))

A similar perspective is to look at the contours of the loss function:

([Source](https://www.quora.com/Why-is-L1-regularization-better-than-L2-regularization-provided-that-all-Norms-are-equivalent))

Optimization Problem¶

Robust PCA can be written:

$$ minimize\; \lVert L \rVert_* + \lambda\lVert S \rVert_1 \\ subject\;to\; L + S = M$$

where:

• $\lVert \cdot \rVert_1$ is the L1 norm. Minimizing the L1 norm results in sparse values. For a matrix, the L1 norm is equal to the maximum absolute column norm.

• $ \lVert \cdot \rVert_* $ is the nuclear norm, which is the L1 norm of the singular values. Trying to minimize this results in sparse singular values --> low rank.

Implementing an algorithm from a paper¶

Source¶

We will use the general primary component pursuit algorithm from this Robust PCA paper (Candes, Li, Ma, Wright), in the specific form of Matrix Completion via the Inexact ALM Method found in section 3.1 of this paper (Lin, Chen, Ma). I also referenced the implemenations found here and here.

The Good Parts¶

Section 1 of Candes, Li, Ma, Wright is nicely written, and section 5 Algorithms is our key interest. You don't need to know the math or understand the proofs to implement an algorithm from a paper. You will need to try different things and comb through resources for useful tidbits. This field has more theoretical researchers and less pragmatic advice. It took months to find what I needed and get this working.

The algorithm shows up on page 29:

needed definitions of $\mathcal{S}$, the Shrinkage operator, and $\mathcal{D}$, the singular value thresholding operator:

Section 3.1 of Chen, Lin, Ma contains a faster vartiation of this:

And Section 4 has some very helpful implementation details on how many singular values to calculate (as well as how to choose the parameter values):

If you want to learn more of the theory:¶

• Convex Optimization by Stephen Boyd (Stanford Prof):
  • OpenEdX Videos
  • Jupyter Notebooks
• Alternating Direction Method of Multipliers (more Stephen Boyd)

Robust PCA (via Primary Component Pursuit)¶

Methods¶

We will use Facebook's Fast Randomized PCA library.

the algorithm again (page 29 of Candes, Li, Ma, and Wright)

Results¶

Extracting a bit of the foreground is easier than identifying the background. To accurately get the background, you need to remove all the foreground, not just parts of it

LU Factorization¶

Both fbpca and our own randomized_range_finder methods used LU factorization, which factors a matrix into the product of a lower triangular matrix and an upper triangular matrix.

Gaussian Elimination¶

This section is based on lectures 20-22 in Trefethen.

If you are unfamiliar with Gaussian elimination or need a refresher, watch this Khan Academy video.

Let's use Gaussian Elimination by hand to review:

A = \begin{pmatrix} 1 & -2 & -2 & -3 \\ 3 & -9 & 0 & -9 \\ -1 & 2 & 4 & 7 \\ -3 & -6 & 26 & 2 \end{pmatrix}

Answer:¶

$$ LU = \begin{bmatrix} 1 & 0 & 0 & 0\\ 3 & 1 & 0 & 0 \\ -1 & 0 & 1 & 0 \\ -3 & 4 & -2 & 1\end{bmatrix} \cdot \begin{bmatrix} 1 & -2 & -2 & -3 \\ 0 & -3 & 6 & 0 \\ 0 & 0 & 2 & 4 \\ 0 & 0 & 0 & 1 \end{bmatrix}$$

Above example is from Lectures 20, 21 of Trefethen.

Gaussian Elimination transform a linear system into an upper triangular one by applying linear transformations on the left. It is triangular triangularization.

$ L_{m-1} \dots L_2 L_1 A = U $

L is unit lower-triangular: all diagonal entries are 1

The LU factorization is useful!

Solving Ax = b becomes LUx = b:

1. find A = LU
2. solve Ly = b
3. solve Ux = y

Work

Work for Gaussian Elimination: $2\cdot\frac{1}{3} n^3$

Memory

Above, we created two new matrices, $L$ and $U$. However, we can store the values of $L$ and $U$ in our matrix A (overwriting the original matrix). Since the diagonal of $L$ is all $1$s, it doesn't need to be stored. Doing factorizations or computations in-place is a common technique in numerical linear algebra to save memory. Note: you wouldn't want to do this if you needed to use your original matrix $A$ again in the future. One of the homework questions is to rewrite the LU method to operate in place.

Consider the matrix $$ A = \begin{bmatrix} 10^{-20} & 1 \\ 1 & 1 \end{bmatrix} $$

By hand, use Gaussian Elimination to calculate what L and U are:

Answer¶

This is the motivation for $LU$ factorization with pivoting.

This also illustrates that LU factorization is stable, but not backward stable. (spoiler alert: even with partial pivoting, LU is "explosively unstable" for certain matrices, yet stable in practice)

Stability¶

An algorithm $\hat{f}$ for a problem $f$ is stable if for each $x$, $$ \frac{\lVert \hat{f}(x) - f(y) \rVert}{ \lVert f(y) \rVert } = \mathcal{O}(\varepsilon_{machine}) $$

for some $y$ with $$ \frac{\lVert y - x \rVert }{\lVert x \rVert} = \mathcal{O}(\varepsilon_{machine}) $$

A stable algorithm gives nearly the right answer to nearly the right question (Trefethen, pg 104)

To translate that:

• right question: $x$
• nearly the right question: $y$
• right answer: $f$
• right answer to nearly the right question: $f(y)$

Backwards Stability¶

Backwards stability is both stronger and simpler than stability.

An algorithm $\hat{f}$ for a problem $f$ is backwards stable if for each $x$, $$ \hat{f}(x) = f(y) $$

for some $y$ with $$ \frac{\lVert y - x \rVert }{\lVert x \rVert} = \mathcal{O}(\varepsilon_{machine}) $$

A backwards stable algorithm gives exactly the right answer to nearly the right question (Trefethen, pg 104)

Translation:

• right question: $x$
• nearly the right question: $y$
• right answer: $f$
• right answer to nearly the right question: $f(y)$

LU factorization with Partial Pivoting¶

Let's now look at the matrix $$ \hat{A} = \begin{bmatrix} 1 & 1 \\ 10^{-20} & 1 \end{bmatrix} $$

By hand, use Gaussian Elimination to calculate what L and U are:

Answer¶

Idea: We can switch the order of the rows around to get more stable answers! This is equivalent to multiplying by a permutation matrix $P$. For instance,

$$\begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} \cdot \begin{bmatrix} 10^{-20} & 1 \\ 1 & 1 \end{bmatrix} = \begin{bmatrix} 1 & 1 \\ 10^{-20} & 1 \end{bmatrix} $$$$ PA = \hat{A} $$

Apply Gaussian elimination to $PA$.

At each step, choose the largest value in column k, and move that row to be row k.

Homework¶

Example¶

Can compare below to answers in Trefethen, page 159:

Partial pivoting permutes the rows. It is such a universal practice, that this is usually what is meant by LU factorization.

Complete pivoting permutes the rows and columns. Complete pivoting is significantly time-consuming and rarely used in practice.

Example¶

Consider the system of equations:

$$ \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 1 \\ -1 & 1 & 0 & 0 & 0 & 1 \\ -1 & -1 & 1 & 0 & 0 & 1 \\ -1 & -1 & -1 & 1 & 0 & 1 \\ -1 & -1 & -1 & -1 & 1 & 1 \\ -1 & -1 & -1 & -1 & -1 & 1 \end{bmatrix} \mathbf{x} = \begin{bmatrix} 1 \\ 1 \\ 1 \\ 1 \\ 2 \\ 1 \end{bmatrix} $$

Exercise¶

Exercise: Let's use Gaussian Elimination on the $5 \times 5$ system.

Scipy has this funtionality as well. Let's look at the solution for the last 5 equations with $n=10,\,20,\,30,\,40,\,50,\,60$.

What is going on when $n=60$?

Theorem: Let the factorization $PA = LU$ of a matrix A be computed by Gaussian Elimination with partial pivoting. The computed (by a computer with Floating Point Arithmetic) matrices $\hat{P}$, $\hat{L}$, and $\hat{U}$ satisfy

$$\hat{L}\hat{U} = \hat{P} A + \delta A, \quad \frac{\delta A}{A} = \mathcal{O}(\rho \varepsilon_{machine}) $$

where $\rho$ is the growth factor,

$$\rho = \frac{max_{i,j} \lvert u_{ij} \rvert }{max_{i,j} \lvert a_{ij} \rvert } $$

For our matrix above, $\rho = 2^{m-1}$

Unstable in theory, stable in practice¶

Stability of most algorithms (such as QR) is straightforward. Not the case for Gaussian Elimination with partial pivoting. Instability in Gaussian elimination (with or without pivoting) arises only if L and/or U is large relative to the size of A.

Trefethen: "Despite examples like (22.4), Gaussian elimination with partial pivoting is utterly stable in practice... In fifty years of computing, no matrix problems that excite an explosive instability are known to have arisen under natural circumstances." [although can easily be constructed as contrived examples]

Although some matrices cause instability, but extraordinarily small proportion of all matrices so "never" arise in practice for statistical reasons. "If you pick a billion matrices at random, you will almost certainly not find one for which Gaussian elimination is unstable."

Further Reading¶

• Gaussian Elimination/LU factorization-- Trefethn Lecture 20
• Pivoting -- Trefethn Lecture 21
• Stability of Gaussian Elimination -- Trefethn Lecture 22

Follow up from last class¶

What is going on with Randomized Projections?¶

We are taking a linear combination (with random weights) of the columns in the matrix below:

It's like a random weighted average. If you take several of these, you end up with columns that are orthonormal to each other.

Johnson-Lindenstrauss Lemma: (from wikipedia) a small set of points in a high-dimensional space can be embedded into a space of much lower dimension in such a way that distances between the points are nearly preserved.

It is desirable to be able to reduce dimensionality of data in a way that preserves relevant structure. The Johnson–Lindenstrauss lemma is a classic result of this type.

History of Gaussian Elimination¶

Fascinating history of Gaussian Elimination. Some highlights:

• First written record of Gaussian elimination from ~200 BC in the Chinese book Nine Chapters on Arithmetic
• The ancient Chinese used colored bamboo rods placed in columns on a "counting board"
• Japanese mathematicican Seki Kowa (16432-1708) carried forward the Chinese elimintion methods and invented the determinant before 1683. Around the same time, Leibniz made similar discoveries independently, but neither Kowa nor Leibniz go credit for their discoveries.
• Gauss referred to the elimination method as being "commonly known" and never claimed to have invented it, although he may have invented the Cholesky decomposition

More history here

Speeding Up Gaussian Elimination¶

Parallelized LU Decomposition LU decomposition can be fully parallelized

Randomized LU Decomposition (2016 article): The randomized LU is fully implemented to run on a standard GPU card without any GPU–CPU data transfer.

Scipy.linalg solve vs lu_solve¶

This problem has a large growth factor $= 2^{59}$. We get the wrong answer using scipy.linalg.lu_solve, but the right answer with scipy.linalg.solve. What is scipy.linalg.solve doing?

Looking at the source code for scipy, we see that it is calling the LAPACK routine gesvx. Here is the Fortran source code for sgesvx (s refers to single, there is also dgesvx for doubles and cgesvx for complex numbers). In the comments, we see that it is computing a reciprocal pivot growth factor, so it is taking into account this growth factor and doing something more complex than plain partial pivot LU factorization.

Block Matrices¶

This is a follow-up to a question about block matrices asked in a previous class. But first,

Ordinary Matrix Multiplication¶

Question: What is the computational complexity (big $\mathcal{O}$) of matrix multiplication for multiplying two $n \times n$ matrices $A \times B = C$?

You can learn (or refresh) about big $\mathcal{O}$ on Codecademy

What this looks like:

for i=1 to n
    {read row i of A into fast memory}
    for j=1 to n
        {read col j of B into fast memory}
        for k=1 to n
            C[i,j] = C[i,j] + A[i,k] x B[k,j]
        {write C[i,j] back to slow memory}

Question: How many reads and writes are made?

Block Matrix Multiplication¶

Divide $A,\, B,\, C$ into $N\times N$ blocks of size $\frac{n}{N} \times \frac{n}{N}$

([Source](http://avishek.net/blog/?p=804))

What this looks like:

for i=1 to N
    for j=1 to N
        for k=1 to N
            {read block (i,k) of A}
            {read block (k,j) of B}
            block (i,j) of C += block of A times block of B
        {write block (i,j) of C back to slow memory}

Question 1: What is the big-$\mathcal{O}$ of this?

Question 2: How many reads and writes are made?

End¶