How to fix theano numpy wrong number of dimensions!

Recently I had to work with Theano, numpy, scipy libraries in python. had to pass hour with this this piece of error which says TypeError: (‘Bad input argument to theano function at index 1(0-based)’, ‘Wrong number of dimensions: expected 2, got 1 with shape (2,).’) So here it goes, how I fixed it.

For the sake helping other fellow programmers, I am posting the whole error log, which will help them via google:

/usr/lib/python2.7/dist-packages/scipy/lib/_util.py:35: DeprecationWarning: Module scipy.linalg.blas.fblas is deprecated, use scipy.linalg.blas instead
DeprecationWarning)
/usr/lib/python2.7/dist-packages/scipy/lib/_util.py:35: DeprecationWarning: Module scipy.linalg.blas.fblas is deprecated, use scipy.linalg.blas instead
DeprecationWarning)
/usr/lib/python2.7/dist-packages/scipy/lib/_util.py:35: DeprecationWarning: Module scipy.linalg.blas.fblas is deprecated, use scipy.linalg.blas instead
DeprecationWarning)
/usr/lib/python2.7/dist-packages/scipy/lib/_util.py:35: DeprecationWarning: Module scipy.linalg.blas.fblas is deprecated, use scipy.linalg.blas instead
DeprecationWarning)
/usr/lib/python2.7/dist-packages/scipy/lib/_util.py:35: DeprecationWarning: Module scipy.linalg.blas.fblas is deprecated, use scipy.linalg.blas instead
DeprecationWarning)
/usr/lib/python2.7/dist-packages/scipy/lib/_util.py:35: DeprecationWarning: Module scipy.linalg.blas.fblas is deprecated, use scipy.linalg.blas instead
DeprecationWarning)
/usr/lib/python2.7/dist-packages/scipy/lib/_util.py:35: DeprecationWarning: Module scipy.linalg.blas.fblas is deprecated, use scipy.linalg.blas instead
DeprecationWarning)
/usr/lib/python2.7/dist-packages/scipy/lib/_util.py:35: DeprecationWarning: Module scipy.linalg.blas.fblas is deprecated, use scipy.linalg.blas instead
DeprecationWarning)
False
Traceback (most recent call last):
File "bgd.py", line 44, in 
theta = gradient_descent_2(alpha, x, y, 2000)
File "bgd.py", line 24, in gradient_descent_2
hypothesis =find_hypothesis(x, theta).T
File "/usr/local/lib/python2.7/dist-packages/Theano-0.6.0-py2.7.egg/theano/compile/function_module.py", line 497, in __call__
allow_downcast=s.allow_downcast)
File "/usr/local/lib/python2.7/dist-packages/Theano-0.6.0-py2.7.egg/theano/tensor/type.py", line 157, in filter
data.shape))
TypeError: ('Bad input argument to theano function at index 1(0-based)', 'Wrong number of dimensions: expected 2, got 1 with shape (2,).')

So basically it is saying that in my code there is a shape mismatch, apparently there should not be any mismatch, because I coded everything, I KNOW what I did! But apparently there is something in numpy that makes a matrix without height, which means with blank height. So it has a shape of (2, ).

Example could be:

import numpy as np
theta = np.array([1,1])
print theta.shape       #(2,)

It works pretty fine when you are working with nupy alone, but it causes tremendous problems when you are using other library like THENO.

Now let me share another snippet, from where I got this error:

from theano import tensor as T, function

T_x,T_y, T_z, T_theta, T_loss = T.dmatrices('x','y','z','theta','loss')
find_hypothesis=function([T_x,T_theta],T.dot(T_x,T_theta))

hypothesis =find_hypothesis(x, theta).T

From code like this you will get the error I previously mentioned.

Now lets talk about how did I fix it, actually it was simple, I passed hours without any reason. You can solve it just by reshapiing your numpy like this:

theta=theta.reshape(theta.shape[0],1)        #had tough time without this,

Maths behind training neural network using Gradient Descent

Well, recently I tried to contribute on a neural network project, I won’t lie I sucked miserably. But it is never too late to learn something and I will be pretty upset if I got another chance to contribute some thing as sexy as neural network and I suck once again. So here I go, learning neural network.

So what does a neural network model look like, it looks like a composition of bunch of layers and every layer contains bunch of neurons. The top most layer is called output layer and the down most layer is called input layer. Depending on the input the data propagates output. There are many many hidden layers resides in between these two. How many hidden layers would you love to have, it is a modeling problem, we need to consider many things to model our neural network. Every neurons are connected with every neurons of its next layer. Every node/neuron has its own activation function, every edge has its own weight, every layer has its own bias. So end of the day, this way or the other every neuron contributes to the output. Now when we are talking about training a neural network we are basically saying, we basically want to set the value of this weight, bias parameter in such a way so that for every input we get the correct result. How do we do that? Thats what this blog is all about.

 

So what can we do? We can go for Empirical risk minimization! So basically we are transforming the training problem to an optimization problem where we need to minimize a loss function of output and desired output. To save papers, or to impress academics we put it that way, argmin \frac{1}{T} \sum l(f(X^{t}, \theta), Y^{t})+\lambda\Omega (\theta), where f() is some function that predicts output, l() is some loss function of predicted output and real outputs. Ω() is a regularization function that takes Θ which is the weights we are predicting. Which has a use to filter on which values we don’t want to take. We will need to smooth our loss function because it is hard to optimize a non-smooth function.

From optimization literature we can use stochastic gradient descent to solve this. Θ= [w_1, b_1,…w_(i+1), b_(i+1)]. So we do what, For N iteration we will find

<img class="mathtex-equation-editor" src="https://web.archive.org/web/20160331013249im_/http://chart.apis.google.com/chart?cht=tx&chl=%5CDelta%20%3D%20-%20%5CDelta_%7B%5CTheta%7D%20l(f(X%5E%7Bt%7D%2C%20%5CTheta%7D)%2C%20Y%5E%7Bt%7D)%20%20-%20%5Clambda%20%5CDelta_%7B%5CTheta%7D%20%5COmega%20(%5CTheta))%0A%5C%5C%5CTheta%20%5Cleftarrow%20%5CTheta%20%2B%20%5Calpha%20%5CDelta" alt="\Delta = – \Delta_{\Theta} l(f(X^{t}, \Theta}), Y^{t}) – \lambda \Delta_{\Theta} \Omega (\Theta))

\\\Theta \leftarrow \Theta + \alpha \Delta” align=”absmiddle” scale=”0″>

This every iteration is known as epoch. So in every epoch we are finding a gradient, now we need to define those functions.

Lets discuss cost function first. As we see our cost function relies on another function f(). F() is basically a prediction function on probablity P(y=c|x), so it gives the probability of y being in class c, when x is given. We would like to maximize this probability, now as we are talking about we are framing this to a minimization problem, so we can define l() as a negative log likelihood problem.l(f(x), y)= - \sum _c 1_{(y=c)} log f(x)_{c} =-log f(x)_{y}
l() is also known as cross entropy of information theory.

Now we will discuss about the gradients, the partial derivatives of our negative log function is
\frac{\delta}{\delta f(x)_{c}} - log f(x)_y =\frac{-1_{(y=c)}}{f(x)_{y}}

As we see, -1 was not a necessary part of the derivatives, but we are adding this as a filter, when y is not c it is 0 so it filters everything else for the term Fc.

So the gradient of the probablity function

<img class="mathtex-equation-editor" src="https://web.archive.org/web/20160331013249im_/http://chart.apis.google.com/chart?cht=tx&chl=%5CDelta%20_%7Bf(x)%7D%20%20-%20log%20f(x)_y%20%3D%20%5Cfrac%7B-1%7D%7Bf(x)_%7By%7D%7D%20%20%5Cleft%20%5B%201_%7B(y%3D0)%7D%20%5C%5C%20.%5C%5C.%5C%5C.%5C%5C%201_%7B(y%3DC-1)%7D%20%20%5Cright%20%5D%5C%5C%20%3D%20%5Cfrac%7B%20-e(y)%20%7D%7Bf(x)_y%7D&quot; alt="\Delta _{f(x)} – log f(x)_y = \frac{-1}{f(x)_{y}} \left [ 1_{(y=0)} \\ .\\.\\.\\ 1_{(y=C-1)} \right ]

\\ = \frac{ -e(y) }{f(x)_y}” align=”absmiddle” scale=”0″>

Now we are interested in partial derivative of output pre-activation function:

\frac{\delta}{\delta a^{(L+1)} X_c} - log f(x)_y = \frac{-1}{f(x)_y}  *  \frac{\delta}{\delta a^{(L+1)} X_c} f(x)_y

Now we replace f() with a softmax function that basically normalizes the exponential of activation it over the summation of other exponentials.

<img class="mathtex-equation-editor" src="https://web.archive.org/web/20160331013249im_/http://chart.apis.google.com/chart?cht=tx&chl=%5Cfrac%7B-1%7D%7Bf(x)_y%7D%20%20*%20%5Cfrac%7B%20%5Cdelta%7D%7B%5Cdelta%20a%5E%7B(L%2B1)%7D%20X_c%7D%20%20softmax%20(a%5E%7B(L%2B1)%7D%20X)_y%5C%5C%5Cfrac%7B-1%7D%7Bf(x)_y%7D%20%20*%20%5Cfrac%7B%20%5Cdelta%7D%7B%5Cdelta%20a%5E%7B(L%2B1)%7D%20X_c%7D%20%20%5Cfrac%7Bexp(a%5E%7B(L%2B1)%7D%20X_y)%20%7D%7B%5Csum%20exp(a%5E%7B(L%2B1)%20%7D%20%20x_c%7D%20%5C%5C&quot; alt="\frac{-1}{f(x)_y} * \frac{ \delta}{\delta a^{(L+1)} X_c} softmax (a^{(L+1)} X)_y\\

\frac{-1}{f(x)_y} * \frac{ \delta}{\delta a^{(L+1)} X_c} \frac{exp(a^{(L+1)} X_y) }{\sum exp(a^{(L+1) } x_c} \\” align=”absmiddle” scale=”0″>

We got this formula for partial derivative of a ratio, \frac{ \delta \frac{g(x)}{h(x)}}{\delta x} } = \frac{ \delta g(x) }{ \delta x } \frac{1}{h(x)} - \frac {g(x)}{h(x)} \frac{\delta h(x)}{\delta x}. If we apply this on our previous equation we get this:

\frac{\delta}{\delta a^{(L+1)} X_c} - log f(x)_y = \frac{-1}{f(x)_y}  *  (1_{(y=c)} softmax(a^{(L+1)}  x_y - softmax(a^{(L+1)}  x_y softmax(a^{(L+1)}  x_c  )

But now we back on our f(x) we got:

\frac{\delta}{\delta a^{(L+1)} X_c} - log f(x)_y = \frac{-1}{f(x)_y}  *  (1_{(y=c) f(x)_y - f(x)_y f(x)_c})

\frac{\delta}{\delta a^{(L+1)} (x)_{c}}} - log f(x)_y =-(1_{(y=c) - f(x)_{c}})

So the gradient,

\Delta  a^{(L+1)}(x) - log f(x)_y = -(e(y)-f(x))–(i)

 

We will also need to find out the gradient of the hidden layers of neural network, if we calculate gradients for each neurons, we will grow old solving this. So we take this equation for chain rule, where a is the activation, p is the loss function, q is the preactivation layer above.

.\frac{\delta}{\delta a} p(a)= \sum_i  \frac{\delta p(a)}{\delta q_i(a)} *  \frac{\delta q_i(a)}{\delta a}

if we are k’th layer we are interested learn their gradients.

Partial derivative at j’th hidden unit and k’th hidden layer, with respect to the activation of hidden layer.

<img class="mathtex-equation-editor" src="https://web.archive.org/web/20160331013249im_/http://chart.apis.google.com/chart?cht=tx&chl=%5Cfrac%7B%5Cdelta%7D%7B%5Cdelta%20h%5E%7Bk%7D%20x_j%7D%20-log%20f(x)_y%3D%20%5Csum_i%20%5Cfrac%7B%20%5Cdelta%20-%20log%20f(x)_y%7D%7B%5Cdelta%20a%5E%7B(k%2B1)%7D%20x_i%20%7D%20%5Cfrac%7B%5Cdelta%20a%5E%7Bk%2B1%7D%20x_i%7D%7B%20%5Cdelta%20h%5Ek%20x_j%20%7D&quot; alt="\frac{\delta}{\delta h^{k} x_j} -log f(x)_y= \sum_i \frac{ \delta – log f(x)_y}{\delta a^{(k+1)} x_i } \frac{\delta a^{k+1} x_i}{ \delta h^k x_j }

” align=”absmiddle” scale=”0″>

a=b+∑W
so differentiation of a will be W, no surprise. Now if we get the gradient of the preactivation it will remain the same.

(W_{,j}^{(k+1)} )^T   \bigtriangledown  a^{k+1} (x) - log f(x)_y

Gradient,

<img class="mathtex-equation-editor" src="https://web.archive.org/web/20160331013249im_/http://chart.apis.google.com/chart?cht=tx&chl=%5Cbigtriangledown%20h%5Ek%20(x)%20-%20log%20f(x)_y%5C%5C(W%5E%7Bk%2B1%7D)%5ET%20(%5Cbigtriangledown_%7B%20a%5E%7Bk%2B1%7D(x)%7D%20-log%20f(x)_y)&quot; alt="\bigtriangledown h^k (x) – log f(x)_y

\\(W^{k+1})^T (\bigtriangledown_{ a^{k+1}(x)} -log f(x)_y)” align=”absmiddle” scale=”0″>—-(2)

<img class="mathtex-equation-editor" src="https://web.archive.org/web/20160331013249im_/http://chart.apis.google.com/chart?cht=tx&chl=%5Cbigtriangledown%20_%20%7Ba%5Ek%20x%20%7D%20-log%20f(x)_y%5C%5C(%5Cbigtriangledown_%7Bh%5Ek%20x%20%7D%20-log%20f(x)_y)%5ET%20(%5Cbigtriangledown_%7Ba%5Ek%20x%20%7D%20h%5Ek%20x)%5C%5C%5C%5C(%5Cbigtriangledown_%7Bh%5Ek%20x%20%7D%20-log%20f(x)_y)%5ET%20%5Cbigodot%20%5B…%2C%20g'(a%5Ek%20x_j)%2C…%5D&quot; alt="\bigtriangledown _ {a^k x } -log f(x)_y

\\(\bigtriangledown_{h^k x } -log f(x)_y)^T (\bigtriangledown_{a^k x } h^k x)

\\\\(\bigtriangledown_{h^k x } -log f(x)_y)^T \bigodot […, g'(a^k x_j),…]” align=”absmiddle” scale=”0″>

Partial derivative of weights:
<img class="mathtex-equation-editor" src="https://web.archive.org/web/20160331013249im_/http://chart.apis.google.com/chart?cht=tx&chl=%5Cfrac%7B%5Cdelta%7D%7B%5Cdelta%20W_%7Bi%2Cj%7D%5Ek%7D%20-log%20f(x)_y%5C%5C%20%3D%5Cfrac%20%7B%20%5Cdelta%20-%20log%20f(x)_y%20%7D%7B%20%5Cdelta%20a%5E%7Bk%7D%20x_i%7D%20*%20%20%5Cfrac%20%7B%20%5Cdelta%20a%5E%7Bk%7D%20x_i%20%7D%7B%20%5Cdelta%20W_%7Bi%2Cj%7D%5Ek%7D&quot; alt="\frac{\delta}{\delta W_{i,j}^k} -log f(x)_y

\\ =\frac { \delta – log f(x)_y }{ \delta a^{k} x_i} * \frac { \delta a^{k} x_i }{ \delta W_{i,j}^k}” align=”absmiddle” scale=”0″>

\\ =\frac { \delta - log f(x)_y }{ \delta a^{k} x_i}  h_j^{k-1} x

which is basically because a (x)=b(k)+∑W(k) h(k-1) (x)

Partial derivatives of Biases:

<img class="mathtex-equation-editor" src="https://web.archive.org/web/20160331013249im_/http://chart.apis.google.com/chart?cht=tx&chl=%5Cfrac%7B%5Cdelta%7D%7B%5Cdelta%20b_%7Bi%7D%5Ek%7D%20-log%20f(x)_y%5C%5C%20%3D%5Cfrac%20%7B%20%5Cdelta%20-%20log%20f(x)_y%20%7D%7B%20%5Cdelta%20a%5E%7Bk%7D%20x_i%7D%20*%20%20%5Cfrac%20%7B%20%5Cdelta%20a%5E%7Bk%7D%20x_i%20%7D%7B%20%5Cdelta%20b_%7Bi%7D%5Ek%7D%5C%5C%3D%5Cfrac%20%7B%20%5Cdelta%20-%20log%20f(x)_y%20%7D%7B%20%5Cdelta%20a%5E%7Bk%7D%20x_i%7D%20*%201&quot; alt="\frac{\delta}{\delta b_{i}^k} -log f(x)_y

\\ =\frac { \delta – log f(x)_y }{ \delta a^{k} x_i} * \frac { \delta a^{k} x_i }{ \delta b_{i}^k}

\\=\frac { \delta – log f(x)_y }{ \delta a^{k} x_i} * 1″ align=”absmiddle” scale=”0″>

We see that somehow most of them has something in common, that is the gradient of pre-activation.

In backward propagation, we assume that we have got f(x) precomputed. So now for all the ks we will find out the gradients for layers. Theta in a 2 hidden layer neural network is defines as [W1,b1,W2,b2] so from my understanding so far, we put it in gradient descent algorithm. So for each iteration we have new value of Ws and bs, which we can put in our gradient descent algorithm.

 

How should we initialize our parameters? Well, looks like it is required a lot of insight to write about it in this blog, but I can say that there is a paper of Glarot and Bengio published in 2010 that suggested to have a gradient H_i,j = Uniform (-b,b). b=�?6/�?(H_k+ H_k-1)

 

Thanks:
1) Martin Thoma, for the inspiration
2) Hugo Larochelle for the tutorial

Notes on Markov Random Field energy function and Reparameterization Example

Markov Random Field is basically a model of joint probability distribution. It is commonly used in many places where an output get biased depending on the output of its neighbors which including Computer Vision. In computer vision one of the major problem that we are aware of is its noise problem. This noise usually depends on its neighbors which satisfies markovian properties.

Now lets represent the neighborhood relationship in a undirected graph where the vertices represents random variables and edges defines the neighborhood. We assume that this random variables that can be picked can be denoted as a finite set L of cardinality h. If I have n random variable then I can say that I have h^n possible combination out of which one is responsible for all. So we can represent it in a probability distribution, p(x). MRF has its own way to assumes this property for p(x). It assumes a vertice conditionally independent of all other vertices other than its neighbor. Mainly because it would help us to use Hammersly Clifford theorem.
MRF1
In graph theory, there is a common term called clique. Clique is basically the vartices that are somehow connected via an edge. If we consider clique potential function of size 2 then clique potential between 2 vartice will return a non negative value. The P(x) is proportional to the multiplication of all clique potentials.
Previously we were busy with only the outputs, but now we will try to predict inputs as well. So x_n were unobserved data and d_n’s are observed data. So now we have another edge which will be subject to clique potential as well. So now the probablity of d given x, P(d|x) is proportional to the multiplication of all clique potentials between x and d.
MRF2
But the function that we get can be normalized using an exponential function, P(x)=Exp(-E(x))/z
Here we get another function, E(x) which is known as Energy function. E(x) is summation instead of product as because we have an exponential in the original function. E(x)=Summation of theta(x)
This E(x) is closely related with Clique Potential because its just analogous to -log(psi(x))
Let me give an example of this energy function. Suppose we are separating out background and foreground image. We can mark which one is backgroun and which one is foreground by marking them in different label. The pixel that has a probability denoted like this P(x_p,l_1)=2 and this P(x_p,l_0)=5 then it is more probable ffor x_p to be in foreground. Now we have to calculate total energy function to get the result. 

E(x)=Summation of theta(x) [for unary potential]
EMFpotential1 

But if we want more accurate result we would go for binary potential.  So we will get connected in all possible way.
E(x)=Summation of theta(xp) + Summation of theta(xp,xq) [for binary potential]

EMFpotential2

Now using dynamic programming we would figure out the minimum energy so that we get a higher probability in total.
Now we know that the only thing we care about is energy function. So we can change the value of potential function to make the calculation more convenient by using Reperamatrization. Reperamatrization is a trick where we change the value of potential function without leting the Energy function being changed. Few of the example of that are given below.

Implementation and mathematics behind Harris Corner Detection

Say I have this amazing friend who travels a lot, say he visits every corner of the world. If I had the job/money I would have done the same, but unfortunately I need to be passionate a little bit to make it happen. So he want to share all his photos with me. There could have been the picture of the hill top buildings of italy, the underground city of turkey, Plitvice lakes of Croatia, Prague of Czech Republic. Obviously some of the photos will be in Zambia capturing The Southern Cross. Bazaruto Archipelago of Mozambique will obviously be there, Some of them will capture Nyiragongo Volcano of Congo. And so on, not one picture in each, but instead, there would be multiples of it. At this age, we don’t just take one picture at a place, right? But as I have seen these places in only in my dream, I could not really tell, which pictures are taken from where just by looking at it. I need to take a closer look, I need to figure out something unique in that picture and correspond that unique feature with other pictures I can tell if the place is same or not. They key point of this para is “Feature detection”.

So as a human we do feature detection. We want our algorithm to do that for us, by nature computer can’t tell a man “a man”, because everything is just a bunch of number that represent pixels, how would it know if they are feature or not? For the simplicity sake lets say we know there are two images where the picture where the image has been captured from the same place. Say we have two corresponding images representing the same scene, we need to match one point of the image to the other point. If we try to correspond each pixel of one image with other image. It will be extremely costly, because even a simple can have 1000×600 pixels. So better approach would be to get few unique identical points that is very rare in the images. Say within 1000×600=100000 pixels, we choose 100 points that does not comes too often, within two images. Since both of the images represent the same scene then something that is unique feature in one image should also be unique in other. So now our new aim will be to correspond 100 points of one image to the other.

Feature pixcels such as mountain peaks, building corners, doorways, or interestingly shaped patches are kinds of localized feature are often called keypoint features, interest points or even corners. They are basically defined by the appearance of patches of pixels surrounding the point location. Another set of important features are edges, say the profile of mountains against the sky.

As we see, similarity plays a great role in it, so we can use summed square difference here between two images, mathematically,

w is the weight, say some pixcel that is close to x_i has a different weight when it is far far away. Usually, it is a gaussian filter.

Which is essentially the correlation equation.

Now we know are going to consider pixels as a series of 2D array, like I(x+u,y+v). As we know any function can be represented as a tylor series, so do our I(). So our I(x+u,y+v) can be represented as

So putting all together:


Now we will compute eigenvalues, it helps us to find least squares plane of some data. It helps us to find the surface normal of point cloud.

is a 2×2 matrix, it can have two eigenvalues, λ1 and λ2.

If λ1 and λ2 are too small then flat.
If λ1 and λ2 are large then it is corner.
If λ1 or λ2 are too small or too large from other it is edge.

Harris measured this cornerness using R=λ1 * λ2 – k (λ1 + λ2)^2 = det M – k(trace M)^2,

if R is negative then it is edge, R is small then flat, R is greater than corner.

from scipy.ndimage import *
from pylab import *
from PIL import *

def compute_harris_response(im,sigma=3):
    """ Compute the Harris corner detector response function for each pixel in a graylevel image. """

    # derivatives
    imx = zeros(im.shape)
    filters.gaussian_filter(im, (sigma,sigma), (0,1), imx)
    imy = zeros(im.shape)
    filters.gaussian_filter(im, (sigma,sigma), (1,0), imy)

    # compute components of the Harris matrix
    Wxx = filters.gaussian_filter(imx*imx,sigma)
    Wxy = filters.gaussian_filter(imx*imy,sigma)
    Wyy = filters.gaussian_filter(imy*imy,sigma)

    # determinant and trace
    Wdet = Wxx*Wyy - Wxy**2
    Wtr = Wxx + Wyy

    return Wdet / Wtr


def get_harris_points(harrisim,min_dist=10,threshold=0.1):
    """ Return corners from a Harris response image
    min_dist is the minimum number of pixels separating
    corners and image boundary. """
    # find top corner candidates above a threshold
    corner_threshold = harrisim.max() * threshold
    harrisim_t = (harrisim > corner_threshold) * 1
    # get coordinates of candidates
    coords = array(harrisim_t.nonzero()).T
    # ...and their values
    candidate_values = [harrisim[c[0],c[1]] for c in coords]
    # sort candidates
    index = argsort(candidate_values)
    # store allowed point locations in array
    allowed_locations = zeros(harrisim.shape)
    allowed_locations[min_dist:-min_dist,min_dist:-min_dist] = 1
    # select the best points taking min_distance into account
    filtered_coords = []
    for i in index:
        if allowed_locations[coords[i,0],coords[i,1]] == 1:
            filtered_coords.append(coords[i])
            allowed_locations[(coords[i,0]-min_dist):(coords[i,0]+min_dist),
            (coords[i,1]-min_dist):(coords[i,1]+min_dist)] = 0
    return filtered_coords

def plot_harris_points(image,filtered_coords):
    """ Plots corners found in image. """
    figure()
    gray()
    imshow(image)
    plot([p[1] for p in filtered_coords],[p[0] for p in filtered_coords],'*')
    axis('off')
    show()






im = array(Image.open('empire.jpg').convert('L'))
harrisim = compute_harris_response(im)

filtered_coords = get_harris_points(harrisim,10)
plot_harris_points(im, filtered_coords)

Thanks to:
1. Professor Mubarak Shah for sharing his course online.
2. Professor Szeleski, for making his book draft public
3. Jan Erik Solem for his book on programming computer vision.
4. That random guy at physics forum who tells me why it works, after full day of wondering, why it works.

Probability transformation, moment and possible uses

This is the second blog on my adventure to learn more about probability and statistics. This is one is about Probability transformation, probability momentum and their use in skewness, kurtosis, variance and so on

Sometime it is easier for us if we can represent a probability distribution in a different way, say we have a very complicated joint probability distribution of a set of random variables, we want to transform it into much simpler form using probability transformation.So, using probability transformation constructs an equivalent set of values that can reasonably be modelled from a specified distribution.

Suppose we have a distribution like this:
F(x) \int_{\sqrt{ \pi /2}}^{0} 2x cos x\^{2} dx

We want to transform it into something else, say y=x^{2}
Replace the limits x = 0 and x = \sqrt{2 / \pi} will change into y=0 and y=\pi/2
\frac{d}{dx}x\^{2}=frac{d}{dx}y =>2x dx=dy
So, F(y) \int_{ \pi /2}^{0} sin y dy will represent the same function.

If cdf F_x(x), and transformation function Y=g(x), then the relationship between F_x(x) and F_y(y) depends on whether g is increasing or decreasing, if it is increasing then F_y(x)=F_x(g\{-1}(y)) but if it is decreasing then F_y(x)=1-F_x(g\{-1}(y))

Let X have pdf fx), let Y=g(x) Suppose there exists a partition A_0, A_1,…A_k of ample saple s[ace
i) g(x)=g_i(x) for x \in A_i
ii)g_i(x) is monotone on A_i
iii) Y=

There are a lot of interesting theorem I found in the book, like say if we have Y=F(x), then P(Y

In binomial transform tation a set of number a_0, a_1, a_2, … get transformed into b_0, b_1, b_2, …, where
b_n =\sum_{k=0}^{n} (-1)^k {(n-k)}\dbinom{n}{k} a_i

A discrete random variable X has a binomial distribution if its pmf is of the form f_x(x)=P(X=x)=\dbinom{n}{k} p^{x} (1-p)^{n-x} where 0< =p=0

If g is increasing function on X, F_y(y)=F_x(g^{-1}(y))
If g is decreasing function on X, F_y(y)=1-F_x(g^{-1}(y))

Now, lets get started with Expectation, what to expect and what not to, this is one hell of a question. Don’t go for very serious issue at the beginning, but at least we can formulate what to expect when we roll a dice. It is calculated as following E(X)=\sum{n}^{n=0} x p(x) So what to expect when we roll the dice? Dice has 6 side, so for each side p(x)=1/6.
E(X)=(1/6)*1+(1/6)*2+(1/6)*3+(1/6)*4+(1/6)*5+(1/6)*6=3.5

So which is basically the mean of the number so the expection is basically the mean, E(X). Not only mean but expectation also has a relation with variance as mean has a relation with variance. Variance \sigma/^{2} can be calculated using the expectation of x-EX )\^2, so E(x-EX )\^2 is variance. So from expectation we know the average, as well as how sparse our data is, on other word variance, so standard deviation is just root of \sigma.

Probability moments are formulated using this eq E(X\^k), so its a expected power of random variable.Central moment E(X -EX)/^{k}. Now we see that anything like mean or variance can be calculated using this probability moment, when k=1 it is same as the equation of mean, when k=2 it is close to variance. For k=3 we E(x-u)\{3}/\sigma\^{3} we get the skew of the moment. skew defines the goodness of the distribution. Skewness=(mean-median)/Standard deviation. Positive distribution pick will be in the left, and vise versa. For k=4 we E(x-u)\{4}/\sigma\^{4} we get the kurtosis, Kurtosis defines whether the distribution has higher peak then the normal distribution or flatter.

If random variable with finite variance, then for any constants a and b, Var(aX+b0=a^2 Var X.

Convergence of MGF theorem states, if lim t->0 M(t)_i=M(t) for all neighbor of 0, then F_xi(x)=F_x(x)

Notes on ensuring Maximum Energy Flow using s-t flow/cut & residual graph

Maximum Flow problem is one of the important problems in combinatorial optimization. It finds feasible solution for flowing energy from one source to another point. Here we are interested to find a minimum energy labeling in a MRF.

Lets define our problem space using directed graph. Every edge denotes a flow and every flow has a capacity (black numbers in diagram) that it can hold. Each time when we make an energy flow we will mark it in the graph in red and obviously it have to be less than the capacity. The excess energy is the difference between incoming and outgoing values. For example E(v1)=1-(4+3)=-6. We can also compute excess energy for a set of veteces, E({v1,v2}). Now we need to consider v1 and v2 as a block and the difference between incoming and outgoing edges must be calculated.  E({v1,v2})=10+1 -(3)=8. We will also get the same result if we compute E(v1) and E(v2) separately and sum them up.flow
Now we will try to explain S-t flow. The main property of a S-t flow is that all the vertice other than source or drain must have an equal value of input and output.
stflow
s-t cut is a way where we separate the graph in two parts in such a order that source and drain goes into the second part and then we find the output capacity difference of it.
stcut
s-t flow can also be represented in term of s-t cut. In terms of s-t flow all other vertices excess energy of 0. So basically the value of flow for s-t flow=E(s) which is similar to write E(s)-E(v) so it is actually an s-t cut.  E(s)-E(v) can’t exceed the capacity.
To solve maximum flow, we can use Residual Graph. Resudual Graph is an undirected graph. In this undirected graph we would remove every edge when it gets saturated in other word we would keep edge of all the graph which has edge less than capacity. So whats our stopping condition? Eventually we will reach to a state where our sink will lose all its edge so its a state where we can’t reach in the end.
rasidual
So we have got the maximum flow for the capacity of that. From the max value min cut theorem we can say that it would give us the minimum cut.

TF-IDF and Cosine similarity on Document Rank Retrieval

I was interested in learning about the algorithms that tries to match document and our queries when it tries to match with the documents and returns on a rank from best match to least match. In this blog I have tried to figure out the question why I have to use tf-idf and cosine similarity .

Jaccard similarity coefficient is too straight forward. For a specific search query It figures out the rank depending on the word it has matched in the total number of unique word it has. I think rather than saying like this, if I say formally then it would be a lot more clear. formally it is defined by the division of cardinality of the set of their intersection and the cardinality of the set of unions.

(1) \begin{equation*} Jaccard(A,B)= \frac{ | A \cup B |}{ | A \cap B |} \end{equation*}

Example:
query=ides o march

document1=caesar died in march
document2=the long march

jaccard(query, document1)=1/6
jaccard(query, document2)=1/5

But Jaccard is too binary to be awesome. In real life if a word frequently appears in a document then the other than the query should retrieval the first document instead of the second one but jaccard coefficient does not count the frequency at all, it checks whether the word appears or not, it is boolean. The problem with Boolean queries are that, it often results too many or too few results. When we use AND it becomes too short and for OR the results become too short. So for long query it becomes too choosy and for short queries it shows everything that matches.
We have detected a problem, so now it is time to solve it. We will try to solve it using Term Frequency Waiting. This time we will try to move away from boolean approach. So this time a presence of a word in a document is not enough but their count is also something of our headache. So the same thing of above but this time we will count the words it has appeared. Suppose for the sentence: “in the month of march, they are going for long march” will vote twice for the document when the query is “march”, in fact it will be 2/11.
tf(query, document)= total query words present in document/ total number of words in query and document

Let me remind the scenario once again so that we can detect any possible bug in our assumption, suppose I am looking for a word and one document has it once and the other document has it 10 times. obviously the 10 time will be the one I would more likely be looking for but it won’t be 10 time more important. So we can use log assumption here. So we can do,

(2) \begin{equation*} w= 1+log(tf) \end{equation*}

But guess what? what is the most used word in our conversation? I am pretty sure that “is, are, it” sort of word will get the most frequency rate but it appears in all the documents which won’t help us to do anything at all. In fact we are interested about rare words – the antique one.
So it was the necessity for introducing Inverse Frequency Model. So let me tell you again the rare terms are more important here for now. For example a document full of the, to is not important to recognize the document anymore because it appears in every document. So in this model the more it appears in the documents, it is less important… in other word the important get divided. So lets do it,

(3) \begin{equation*} idf=\frac{N}{df} \end{equation*}

where N= total number of documents we have,
df= the number of document that has our query word
But again it does not increase like a plain line, so lets use logarithmic scale for the same reason as above.

(4) \begin{equation*} idf=Log \frac{N}{df} \end{equation*}

more formally:

(5) \begin{equation*} idf(t,D)= log(tf) \end{equation*}

so now for 1 occurance in 1M will get idf = 6

TF only care about the number of time it appears in the document and idf cares about the rare words. But we need a proper balance of things. It should be rare but at the same time that rare word should appear multiple time in the document. So how would we balance? we can just multiply both of them. tf-idf is the best weighting scheme known.
Let me help you a little bit more with its implementation. When you are done your tf-idf would look like a 1xm vector matrix. somewhat like array.
[It, is, March]=[0.33,0.33,0]
Remember how do we draw vectors? if a vector is represented as [x,y] just ploting line that go through [0,0] and that connects [x,y]. TF-IDF makes it multidimensional nothing else. we have represented something using tf-idf, thats cool, but how would we we retrive ranks? We can make tf-idf matrix for both query and document then calculate vector distance. But here we can face a new problem because of its vector nature. remember that [2,2] and [4,4] represent the same line, just twice as large as the line. Which is a drawback of this tf-idf because it adds the distance between query and document very high. So to get out of this draw back we will introduce cosine similarity. Basically cosine similarity is pretty simple idea, instead of distance now we will calculate the angle between the lines.

(6) \begin{equation*} cos(\theta)=\frac{A \dot B}{| A | | B |} \end{equation*}

so from angle now we can detect the similarity.