Randomized SVD demos¶

Part 1: Random Projections (with word vectors)¶

The purpose of this section is to illustrate the idea of random projections preserving structure with the concrete example of word vectors!

To use language in machine learning (for instance, how Skype translator translates between languages, or how Gmail Smart Reply automatically suggests possible responses for your emails), we need to represent words as vectors.

We can represent words as 100 dimensional vectors, using Google's Word2Vec or Stanford's GloVe. For example, here is the word "python" as a vector in GloVe:

Goal: Use randomness to reduce this from 100 dimensions to 20. Check that similar words are still grouped together.

More info: If you are interested in word embeddings and want more detail, I gave a longer workshop about them available here (with code demo).

Style note: I use collapsible headings and jupyter themes

Loading the data¶

The dataset is available at http://files.fast.ai/models/glove/6B.100d.tgz To download and unzip the files from the command line, you can run:

wget http://files.fast.ai/models/glove_50_glove_100.tgz 
tar xvzf glove_50_glove_100.tgz

You will need to update the path below to be accurate for where you are storing the data.

What our data looks like¶

We have a long list of words:

wordidx allows us to look up a word in order to find out it's index:

Words as vectors¶

The word "python" is represented by the 100 dimensional vector:

This lets us do some useful calculations. For instance, we can see how far apart two words are using a distance metric:

Smaller numbers mean two words are closer together, larger numbers mean they are further apart.

The distance between similar words is low:

And the distance between unrelated words is high:

Bias¶

There is a lot of opportunity for bias:

I just checked the distance between pairs of words, because this is a quick and simple way to illustrate the concept. It is also a very noisy approach, and researchers approach this problem in more systematic ways.

I talk about bias in much greater depth in this workshop.

Visualizations¶

Let's visualize some words!

We will use Plotly, a Python library to make interactive graphs (note: everything below is done without creating an account, with the free, offline version of Plotly).

Methods¶

This method will pick out the 3 dimensions that best separate our categories from one another (stored in dist_btwn_cats), while minimizing the distance of the words within a given category (stored in dist_within_cats).

Preparing the Data¶

Let's plot words from a few different categories:

Again, we need to look up the indices of our words using the wordidx dictionary:

Now, we will make a set combining our words with the first 10,000 words in our entire set of words (some of the words will already be in there), and create a matrix of their embeddings.

Viewing the words in 3D¶

The words are in 100 dimensions and we need a way to visualize them in 3D.

We will use Principal Component Analysis (PCA), a widely used technique with many applications, including visualizing high-dimensional data sets in a lower dimension!

PCA¶

Random Projections¶

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 (proof uses random projections).

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

Part 2: Background Removal¶

Our goal today:

Load and Format the Data¶

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

Import needed libraries:

Helper Methods¶

Format the Data¶

An image from 1 moment in time is 120 pixels by 160 pixels (when scaled). We can unroll that picture into a single tall column. So instead of having a 2D picture that is $120 \times 160$, we have a $1 \times 19,200$ 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 hundredth of a second for 100 seconds (so 10,000 different images, each from a different point in time), we'd have a $10,000 \times 19,200$ 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.

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

Part 2 continued: SVD¶

Intro to SVD¶

Applications of SVD:

• Principal Component Analysis
• Data compression
• Pseudo-inverse
• Collaborative Filtering
• Topic Modeling
• Background Removal
• Removing Corrupted Data

“a convenient way for breaking a matrix into simpler, meaningful pieces we care about” – David Austin

“the most important linear algebra concept I don’t remember learning” - Daniel Lemire

This is really slow, so you may want to save your result to use in the future.

In the future, you can just load what you've saved:

What do $U$, $S$, and $V$ look like?

Check that they are a decomposition of M

They are! :-)

s is the diagonal of a diagonal matrix

Do you see anything about the order for $s$?

$U$ is a giant matrix, so let's just look at a tiny bit of it:

Finding the background¶

How do we get the people from here?

High-resolution version

Make Video¶

I was inspired by the work of fast.ai student Samir Moussa to make videos of the people.

Speed of SVD for different size matrices¶

2 Backgrounds in 1 Video¶

We'll now use real video 008 dataset from BMC 2012 Background Models Challenge Dataset, in addition to video 003 that we used above.

This is slow, so you may want to save your result to use in the future.

In the future, you can just load what you've saved:

The top few components of U:

Background removal:

Aside about data compression¶

Suppose we take 700 singular values (remember, there are 10000 singular values total)

space saved = data in U, s, V for 700 singular values / original matrix

We only need to store 15.3% as much data and can keep the accuracy to 1e-5! That's great!

Now back to our background removal problem:¶

That's pretty neat!!! but...¶

Downside: this was really slow (also, we threw away a lot of our calculation)

The runtime complexity for SVD is $\mathcal{O}(\text{min}(m^2 n,\; m n^2))$

Let's create a smaller version of M¶

Idea: Let's use a smaller matrix!

We haven't found a better general SVD method, we'll just use the method we have on a smaller matrix.

How do we get the people from here?

What this method is doing¶

And the people:

Timing Comparison¶

Full SVD:

Our (overly simplified) randomized_svd from above:

Scikit learn:

Randomized SVD from Facebook's library fbpca:

I would choose fbpca, since it's faster than sklearn but more robust and more accurate than our simple implementation.

Here are some results from Facebook Research:

Time and Accuracy vary with k¶

Math Details¶

Process behind Randomized SVD¶

Here is a process to calculate a truncated SVD, described in Finding Structure with Randomness: Probabilistic Algorithms for Constructing Approximate Matrix Decompositions and summarized in this blog post:

1. Compute an approximation to the range of $A$. That is, we want $Q$ with $r$ orthonormal columns such that $$A \approx QQ^TA$$

2. Construct $B = Q^T A$, which is small ($r\times n$)

3. Compute the SVD of $B$ by standard methods (fast since $B$ is smaller than $A$), $B = S\,\Sigma V^T$

4. Since $$ A \approx Q Q^T A = Q (S\,\Sigma V^T)$$ if we set $U = QS$, then we have a low rank approximation $A \approx U \Sigma V^T$.

So how do we find $Q$ (in step 1)?¶

To estimate the range of $A$, we can just take a bunch of random vectors $w_i$, evaluate the subspace formed by $Aw_i$. We can form a matrix $W$ with the $w_i$ as it's columns. Now, we take the QR decomposition of $AW = QR$, then the columns of $Q$ form an orthonormal basis for $AW$, which is the range of $A$.

Since the matrix $AW$ of the product has far more rows than columns and therefore, approximately, orthonormal columns. This is simple probability - with lots of rows, and few columns, it's unlikely that the columns are linearly dependent.

Why M ~ Q Q.T M¶

We are trying to find a matrix Q such that $M \approx Q Q^T M$. We are interested in the range of $M$, let's call this $MX$. $Q$ has orthonormal columns so $Q^TQ = I$ (but $QQ^T$ isn't $I$, since $Q$ is rectangular)

so... $$ MX = QQ^TMX $$

If $X$ is the identity, we'd be done (but then $X$ would be too big, and we wouldn't get the speed up we're looking for). In our problem, $X$ is just a small random matrix. The Johnson-Lindenstrauss Lemma provides some justification of why this works.

The QR Decomposition¶

We will be learning about the QR decomposition in depth later on. For now, you just need to know that $A = QR$, where $Q$ consists of orthonormal columns, and $R$ is upper triangular. Trefethen says that the QR decomposition is the most important idea in numerical linear algebra! We will definitely be returning to it.

How should we choose $r$?¶

Suppose our matrix has 100 columns, and we want 5 columns in U and V. To be safe, we should project our matrix onto an orthogonal basis with a few more rows and columns than 5 (let's use 15). At the end, we will just grab the first 5 columns of U and V

So even although our projection was only approximate, by making it a bit bigger than we need, we can make up for the loss of accuracy (since we're only taking a subset later).

How this is different from random mean¶

Random mean:

Mean image:

End¶