Bishop PRML Chapter 1 (1.15-1.16)

Finding \( n(D,M) \) (problem 1.15) is finding the number of ways we can sum \(D\) non-negative numbers to M. This is a basic combinatorics problem with a number of interesting solution methods. I will outline three approaches that feel somewhat different and find application in solving other problems.

\[ \mathrm{Solve:} \quad x_1+x_2+...+x_D=M \]
\[ x_i \gt =0 \quad 1 \le i \le D\]

Method 1: Using generating functions

Coefficient of \(x^M\) in \( \left(1+x+x^2+....\right)^D \)

=Coefficient of \(x^M\) in \( \left(1-x\right)^{-D} \)

Using -ve binomial  series the coefficient is \( \left( \begin{matrix} D+M-1\\M  \end{matrix} \right) \)

Method 2: Bars and Balls

Lets say we have a set of \(M\) 1’s, and a set of \(D-1\) bars. The \(D-1\) bars split the 1’s into the number of 1's representing each of the \(D x_i\)’s . So in total we have \(D+M-1\) positions out of which we select \(M\) position for 1’s and the rest are for the bars.

A following similar approach yields the wrong solution: So we have \(M\) 1’s so we can select each bar in \(M+1\) ways. The mistake here is that selecting the first and second bars in positions 1 and 2 is the same as selecting them in positions 2 and 1. We need each new bar to come to the right of the ones we have already positioned. Okay, so can't we just select the position of \(D-1\) bars from \(M+1\) positions. This gives us the solution to the problem where each \(x_i \gt 0\) because we cannot choose the same position for two bars. However, when \(x_i \ge 0\) we are allowed to choose bars from the same position between balls so we start off with a set of positions for bars and balls and then select the positions of the bars.

Method 3: Indices

Let each \(x_i\) represent the number of balls in box \(i\). The number of boxes is \(D\) and the number of balls is \(M\). Represent each solution as a non-decreasing list of box numbers. E.g. 11335 represents 2 balls in box 1, two balls in box 3 and one ball in box 5. Now create another sequence which is the indices of the balls 01234. Add the two sequences 12569. We can see that all increasing sequences can be generated with a smallest value of 1 and a largest value of \(D+M-1\). So the number of ways is again \( \left( \begin{matrix} D+M-1\\M  \end{matrix} \right) \).

The extension to problem 1.16 is then fairly straightforward. We just add another variable giving us a total of \(D+1\) dimensions and solve the same problem as above. The newly added variable will take over the leftover whenever the other \(D\) variables sum to less than \(M\). So, in this case we have the number of ways as \( \left( \begin{matrix} D+M\\M  \end{matrix} \right) \)

Puzzles: Languages and locks

I recently heard a fun puzzle. Puzzle: How many languages are needed so that we can have that for any two people in a group of 70 each one knows a language that the other does not.






Solution: Basically we need to assign to each person a unique subset of languages. If there are K languages number of such such subsets will be maximum if each person knows K/2 languages. So, we will need the number of languages K to be such that

\[ \left( \begin{matrix} K\\K/2  \end{matrix} \right) \ge 70 \]

Another, variant of this puzzle is regarding bank locks. Puzzle: Suppose, you have \(n\) people working at a bank and you want that at least \(m\) should be present to be able to open a lock. How many locks do you need, and how many keys are needed for each lock.




Solution: For each subset of \(m-1\) people there should be a lock that they cannot open so  we need

\[ \left( \begin{matrix} n\\m-1  \end{matrix} \right) \]

locks. Also, for any lock we cannot find m people that together cannot open it. So each lock should have \(n-(m-1)\) keys.

Trace and transpose trick: Computing derivatives

I am not sure why but none of the college courses on linear algebra or calculus ever introduced me to working with derivatives of functions of matrices and vectors. This comes up on a daily basis in optimization and is a lot cleaner than dealing with summations. Here is a trick that comes in very handy when taking derivatives of a functions that maps vectors or matrices to scalars. Basically, one can replace a scalar results of matrix vector products by its transpose and in some cases by their trace (and then use the cyclic property of  trace to rewrite expressions) in which the desired terms appear at the end and are correctly transposed.

Simplest case: Linear function takes a vector to a scalar.
\[f : \mathbb{R}^n \rightarrow \mathbb{R}\]
\[ \mathrm{d} \left(c^Tx \right) = c^T \mathrm{d}x\]

Next level: Function takes a vector to a scalar using a quadratic form. We will use the transpose trick here:
\[f : \mathbb{R}^n \rightarrow \mathbb{R}\]
\[ \mathrm{d} \left(x^TAx \right) = (\mathrm{d}x)^TAx + x^TA(\mathrm{d}x) \]

The first term on the RHS is a scalar so we can take its transpose and get

\[ \mathrm{d} \left(x^TAx \right) = x^TA^T(\mathrm{d}x) + x^TA(\mathrm{d}x)  = (x^TA^T + x^TA) \mathrm{d}x\]

Next level: Function takes a matrix to a scalar. We will use the trace trick here: 

\[f : \mathbb{R}^{nxn} \rightarrow \mathbb{R}\]
\[ \mathrm{d} \left(a^TXb \right) = a^T (\mathrm{d}X)b \]
Now since the RHS is a scalar we can think of it as a trace and rotate the terms to get:
\[ \mathrm{d} \left(a^TXb \right) = \mathrm{Tr}(ba^T \mathrm{d}X) \].
Now, we will use the symbol (:) that rewrites the matrix as a concatenation of its vectors and rewrite the above as:
\[ \mathrm{d} \left(a^TXb \right) = (ab^T): (\mathrm{d}X): \].
Giving us \( \frac{\mathrm{d} \left(a^TXb \right)}{\mathrm{d}X} = ab^T \).

Lets use this machinery to attempt something more complicated and find \( \frac{\mathrm{d}(a^TX^TCXb)}{\mathrm{d}X} \)

\[ \mathrm{d} (a^TX^TCXb)  = a^T(\mathrm{d}X^T)CXb + a^TX^TC(\mathrm{d}X)b \]

\[ = b^TX^TC^T(dX)a + a^TX^TC(dX)b \]

\[ = \mathrm{Tr}( b^TX^TC^T(\mathrm{d}X)a + a^TX^TC(\mathrm{d}X)b ) \]

\[ = \mathrm{Tr}( ab^TX^TC^T(\mathrm{d}X) + ba^TX^TC(\mathrm{d}X)) \]

\[ = (CXba^T+C^TXab^T): (\mathrm{d}X): \]

\[ \Rightarrow \mathrm{d}(a^TX^TCXb)/\mathrm{d}X = (CXba^T+C^TXab^T) \]

Bell states to test nature's inherent uncertainity

This post is to succinctly describe the key ideas that are presented by Prof. Vazirani in his QMQC course on EdX (I saw this when it was originally offered on coursera) related to Bell's idea to show that quantum randomness is (or is not) inherent in nature.

So, as usual, we have two players Alice and Bob who can initially agree on a strategy and are then separated by a huge distance over which we will assume that they cannot communicate. Now, the game is that both Alice and Bob receive a binary input and produce a binary input. They win when they generate different inputs only when both their inputs are 1, otherwise they should produce the same output. Classically, we can see that the best that Alice and Bob can do is win with a probability of 0.75.

However, if they can share an entangled Bell state then they can achieve a better probability of success. An entangled state of two particles is one which cannot be expressed as a product of the states of the two particles. So, as an example,

\[ \lvert \Psi \rangle = \frac{1}{\sqrt{2}} \left(   \lvert 00 \rangle +  \lvert 11 \rangle \right) \]

is an entangled state. It is instructive to see that this state cannot be written as a tensor product of the individual states of the two particles. Now, the next fun thing to note is that the representation of this state when written in any other orthonormal basis \(u, u^{\perp} \) does not change (again, expand out the tensor products to see this). So, our state above can be expressed as


\[ \lvert \Psi \rangle = \frac{1}{\sqrt{2}} \left(   \lvert uu \rangle +  \lvert u^{\perp}u^{\perp} \rangle \right) \]

Now, let Alice and Bob share the entangled state with one particle each before separation and are then taken to two corners of the universe. They can choose their basis (x for Alice and y for Bob) as shown in the figure below and announce their result depending on whether they see the particle oriented along their basis vector or orthogonal to it. In, all cases except when \(x=y=1\) their basis vectors are oriented at an angle of \( \frac{\pi}{8} \), and only when \(x=y=1\) there basis vectors are orientated at an angle of \( \frac{3\pi}{8} \). In both cases the probability of winning is \( \cos^2\frac{\pi}{8} \approx 0.85 \), much better than can be obtained classically.

A nice article discussing the questions still surrounding the experimental results with Bell's states appeared here recently.

Linear Regression

I recently gave a talk on Linear Regression and its extensions (Ridge regression, LASSO, LARS etc.) as part of the Austin SIGKDD Advanced Machine Learning Meetup. The slides can be found here.

One of the interesting things that I learnt about while preparing for this talk was the fact that coordinate descent works for function where the non-differentiable part of the objective function is separable (slides 29-30). This is why we can use gradient optimization for minimizing the LASSO objective function. The proof works out and the only background that is needed is subgradients. Also, it is instructive to make sure that the proof "really" uses the property of separability and why the proof breaks when we do not have separability.

Hint: For a separable function the subgradient set is a direct product of the subgradient sets of the individual components. However, if the non-differentiable function is not separable then the negative gradient of the differentiable part may not lie in the subgradient set of the non-differentiable function.

Another, interesting thing was the connection between LASSO and LARS. LASSO minimizes the quadratic residual norm with a \(\mathrm{L1}\) penalty and requires that the sign of the coefficient related to a vector be the same as the sign of the inner product of the residual with the vector (easy to show using the Lagrangian). LARS does not enforce this and can therefore produce different solutions. LARS can modified to drop a vector from the active set when this sign change happens and can thus be tailored to generate a LASSO solution. However, then the number of steps taken by LARS is not bounded by \(n\), the dimensionality of the problem.

PCA vs Least Squares

We all know that least squares minimizes the sum of the squared errors. The error is defined as the error between the measured and the predicted value of the dependent variable, where the predicted value is a linear function of the independent variables. If both the dependent and the independent variable are measured quantities and may be noisy, then we might be more interested in minimizing the error of the orthogonal projection of the data points on a linear hyperplane, and this turns out to be the first principal component of the data. The page here has a couple of good images that get this point across. But, lets see that is indeed the case, again, with matrix algebra :)

Let's start with the basics. The projection of a vector \(x\) on a vector \(w\) can be expressed as \(w^Tx \frac{w}{\| w \|^2}\), which in the case of a unit-norm vector \(w\) reduces to \((w^Tx) w\). So, if we have a data matrix \(X\), where the columns of the data matrix form the data points, then the projection of the entire data matrix can be written as \((w^TX) w\). Basically, the projections are all pointing along \(w\) with lengths given by \(w^TX\). In PCA, we seek to find a direction that has the maximum spread (or variance) and we can continue to consider the direction as a vector through the origin if we assume that the data matrix \(X\) has \(0\) mean. Equivalently, we want to find a direction \(w\) that maximizes the following:

\[\mathrm{Cov}\left( w^TX\right) = w^T\mathrm{Cov}\left(X\right)w = w^TXX^Tw\]

Now, lets see what we get for the sum of the squared orthogonal error of the data points on to the line given by the vector \(w\). Each data points has an error given by the 2-norm of the vector that represents the orthogonal error. The sum of the squared errors for each data point can then be represented by a Frobenius norm of a matrix that is composed of these error vectors.

\[\|X - w (w^TX)\|_F = \mathrm{Tr}( (X-ww^TX)^T (X - ww^TX) )\]

\[= \mathrm{Tr}\left(X^TX - X^Tww^TX-X^Tww^TX+X^Tww^Tww^TX\right)\]

\[= \mathrm{Tr}\left(X^TX - X^Tww^TX-X^Tww^TX+X^Tww^TX\right)\] (as \(w\) is unit norm)

\[=\mathrm{Tr}\left( X^TX \right) - \mathrm{Tr}\left(X^Tww^TX\right)\]

\[=\mathrm{Tr}\left( X^TX \right) - \mathrm{Tr}\left(w^TXX^Tw\right)\] (Since trace of cyclic rotations of matrix products is the same)

\[=\mathrm{Tr}\left( X^TX \right) - \left(w^TXX^Tw\right)\] (The trace of 1x1 matrix is the matrix itself)

which is the same as the case above for PCA except for a constant term (\(\mathrm{Tr} (X^TX)\) that is independent of \(w\)) and that we seek to maximize the PCA variance and minimize the orthogonal error. Now, the matrix \(A=XX^T\) represents the covariance matrix of the data points and one seeks to find a unit-norm direction that minimizes the quadratic form \(w^TAw\) for a positive semi-definite symmetric matrix \(A\). One can either construct the Lagrangian or write \(w\) in terms of the eigenbasis to see that we need to choose \(w\) as the eigenvector with the largest eigenvalue.

Logistic Regression: Gradients computed the matrix way

I really enjoy performing calculus with matrix/vector operations rather than messing around with tons of summations. The other nice thing about calculations with matrices/vectors is that the results come out cleanly in a form that can be used directly in Matlab or a similar high-level language and where the matrix/vector multiplication based implementations are typically much faster than using your own for-loops to compute the result. One of the nice examples of this approach is computing the gradients for logistic regression without using any summations at all!

Before we get to logistic regression let's build up one small piece of the machinery (vectorized functions) that will help us when we get there. Lets begin with vectors \(p, z \in \mathbb{R}^n\) and let \(f_v\) be a vectorized version of function \(f\) that applies the function \(f\) component-wise.

\[p = f_v(z)\]
\[\mathrm{d}p = \mathrm{Diag}((f')_v(z))\mathrm{d}z\].

Now, our logistic regression objective function is

\[J(\theta) = \frac{1}{m} \left(   -y^T(log(p)) -(1-y)^T(log(1-p))   \right)\]

where

\[p=g_v(z)\]
\[g(x)=\frac{1}{1+\exp^{-x}}\]
\[z=X^T\theta\]

Lets start with the simplest derivatives and then use those to compute the gradient of the objective function

\[\mathrm{d}g=\frac{\exp^{-x}}{(1+\exp^{-x})^2}\mathrm{d}x=g(x)(1-g(x))\mathrm{d}x\]

Now, lets get to \(p\) and make use of machinery we developed at the very beginning:

\[\mathrm{d}p = \mathrm{Diag}((g')_v(z)) \mathrm{d}z\]
\[\mathrm{d}p = \mathrm{Diag}(g_v(z)(1-g_v(z))\mathrm{d}z=\mathrm{Diag}(g_v(z)(1-g_v(z))X^T\mathrm{d}\theta\]

Now, we come to our original goal of computing the gradient of the objective function \(J\)

\[\mathrm{d}J = \frac{1}{m} \left(   -y^T(\mathrm{Diag}(1/p)) + (1-y)^T(\mathrm{Diag}(1/(1-p)))   \right)\mathrm{d}p\]

\[\mathrm{d}J = \frac{1}{m} \left(   -y^T(\mathrm{Diag}(1/g_v(z))) + (1-y)^T(\mathrm{Diag}(1/(1-g_v(z))))   \right)\mathrm{Diag}(g_v(z)(1-g_v(z))X^T\mathrm{d}\theta\]

\[\mathrm{d}J = \frac{1}{m} \left(   -y^T(\mathrm{Diag}(1-g_v(z))) + (1-y)^T(\mathrm{Diag}(g_v(z)))   \right)X^T\mathrm{d}\theta\]

\[\mathrm{d}J = \frac{1}{m} \left(   g_v(z)^T -y ^T    \right)X^T\mathrm{d}\theta\]

So, finally we have

\[\nabla J = X(g_v(z)-y)\]

A good reference for matrix calculus can be found here.