Showing posts with label matrix algebra. Show all posts
Showing posts with label matrix algebra. Show all posts

Monday, February 10, 2014

Machine learning in a few words

Machine learning is becoming a buzzword, everybody talks aboit it and few seem to be interested in the math underneath (I find statements like "I wanted to know more but all sources were too statistical/mathematical and I wanted more practical stuff").

Let me tell you something:
First: You can't really use Machine Learning if you don't know the statistical/mathematical basis
Second: You can't really use Machine Learning if you don't know the statistical/mathematical basis
Third: You can't really use Machine Learning if you don't know the statistical/mathematical basis

Machine Learning is just a fancy word for the statistical/mathematical tools lying underneath, whose objective is to extract something that we may loosely call knowledge (or something that we understand) from data (or something chaotic that we do not understand), so that computers may take action based on the inferred knowledge. An example of this would be a robot arm/humanoid: without programming actions on direction/velocity/acceleration vectors based on an established model, we may put sensonrs on a subject's articulations, and from these datapoints learn a regression model on the manifold of natural movements. Another example is in Business Intelligente: We may learn groups of customers (market segmentation) so that we may engage several groups with specific policies or offers target at them.

Maching Learning is applied Statistics/Mathematics. Is very little and very unpractical without Optimization/Operations Research, from the algorithmic and practical/scalable point of view.

I've come to the conclusion that there exists two large main approaches to ML, disregarding the specific technique we are dealing with and its target (i.e., supervised or unsupervised), plus one in the middleway:
  • Functional approach (Mathematical)
  • Neural Network/Deep Learning approach (Middle way)
  • Probabilistic approach (Statistical)
In the functional approach, one uses the theory of Hilbert spaces (and therefore of differential equations and linear algebra). The goal is to find a set of transformations of the data so as to best perform in a score for the task (called functional). These transformations come from a pre-defined set which is not related to the data in any way and are a combination of possibly orthogonal basis of a space of functions (transformation) defined on the domain of the original data. Examples of this are: Linear/Ridge/Sparse Regression (linear or identity transformation for regression), SVM (non-linear using the kernel trick), PCA (SVD/eig computation) and KPCA, Matrix factorizations (for signal separation/clusering...), K-means, Projection pursuit... The basic idea is:
I have my data set and a buch of (linear or non-linear) transformations, find a solution applying these transformations to my data so that I can maximize a score functional that I like for my problem. If I want to predict something (classification/regression) then combine these transformations so that the combination fits best my target variable. If I want to examine the nature of the data (dimensionality reduction, matrix factorization, clustering), then use a combination of my transformations to lose as least information as possible (measured with a functional, again).
In the probabilistic approach, the prior knowledge is an assupmtion of the prior probability and the likelyhood, and works towards obtaining a posterior probability (that the outcome is a given choice given the data just seen). Examples are: Logistic Regression (simple non-linear transformation for classification), Naive Bayes (classification), SNE, Gaussian Processes... The general idea is
I want to see independence at the end of the process, therefore I can assume Gaussian multivariate variables so that a linear or non-linear transformation gives me components that are as independent as possible (for unsupervised learning) or assume a probability distribution at the output and a likelyhood and compute a model so that I best fit the likelihood of seeing the target variable given the data.
Neural networks and deep learning are a different story. I consider them to be their own field, drawing tools from the above. They are neither functional because they are not dealing directly with functions (transformations) in the functional analytic setting, and they are not probabilistic for obvious reasons, but use any probability or information theoretic tool as needed. The fact that they connect output of transformations to inputs makes this field related to the first approach indeed, seen as a chain of transformations (spaces defined on the image of their predecessors), but the focus here is obviously the algorithms and that changes things.

Reinforment learning is nowadays just a fancy word for techniques widely known, studied and used in Stochastic Processes, HMM (Hidden Markov Models) being the only exception. Sometimes they call it Sequential Learning, but it is not widely considered Machine Learning, neither scholarly nor popularly.

So, as you can see, there is nothing new (not at least as the discovery of the fundamental theorem of calculus, or of quantum mechanics).

Regarding novelty, I am annoyed each time I read about ML techniques from the big data guys, who are normally programmers starting to do something more that data aggregation and querying. It seems that there is something both new and exotic in what they are saying, and it it mostly well known techniques from statistics, such as this article, or the beginning of Sean Owen's presentation at Big Data Spain.

Now, the practical side of things require that ML scales to Big Data. That limits the applicability of matrices and non-linear transformations to adapt linear methods, so let's see what the next breakthrough is.

Saturday, February 9, 2013

Neil Lawrence's opening course

Neil Lawrence is Professor at the University of Sheffield. He has worked on unsupervised learning for a long time, and has developed algorithms applicable to dimensionality reduction such as the Gaussian Process Latent Variable model (GP-LVM), see the JMLR paper here.

He has a superb inaugural lecture in which he talks about Machine Learning. The link to the starting page to see that video is here. It opens some embedded and annoying player but it is worth dealing with it.

Neil Lawrence's Inaugural Lecture

Title: Life, The Universe and Machine Learning

Time: 17:15 Thursday 6th September 2012

Venue: St George's Church Lecture Theatre, University of Sheffield

Abstract
What is Machine Learning? Why is it useful for us? Machine learning algorithms are the engines that are driving forward an intelligent internet. They are allowing us to uncover the causes of cancer and helping us understand the way the universe is put together. They are suggesting who your friends are on facebook, enabling driverless cars and causing flagging potentially fraudulent transactions on your credit card. To put it simply, machine learning is about understanding data.
In this lecture I will try and give a sense of the challenges we face in machine learning, with a particular focus on those that have inspired my research. We will look at applications of data modelling from the early 19th century to the present, and see how they relate to modern machine learning. There will be a particular focus on dealing with uncertainty: something humans are good at, but an area where computers have typically struggled. We will emphasize the role of uncertainty in data modelling and hope to persuade the audience that correct handling of uncertainty may be one of the keys to intelligent systems.

Monday, October 1, 2012

Jacobi eigenvalue method implementation in C++


The Jacobi eigenvalue problem is an algorithm to compute the eigenvalues of a matrix by canceling out the off-diagonal elements by multiplying the matrices with rotation matrices.

I started with the code in Wikipedia, which is essentially wrong.

Here I offer an implementation that uses Boost uBLAS matrices. It uses a trick to avoid computing sines and cosines. It needs to be improved to avoid multiplication of full matrices when doing the rotations.

This code is copylefted, use it as you wish.

#include "stdafx.h"

#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>

using namespace boost::numeric::ublas;
using namespace boost;

#define EPS 10e-08


void abs(matrix<double> &M)
{
 int n = M.size1();

 for (int k=0; k<n; k++)
 {
  for (int i=0; i<n;i++)
  {
   M(k,i)=abs(M(k,i));
  }
 }
}


void findMax(matrix<double> &M, int &row, int &col)
{
 double m = M(0,0);
 row = col = 0;
 int n = M.size1();

 for (int k=0; k<n; k++)
 {
  for (int l=0; l<n;l++)
  {
   if(M(k,l) > m)
   {
    row = k;
    col = l;
    m = M(k,l);
   }
  }
 }
}


void zeroRowCol(matrix<double> &S, matrix<double> &U, int row, int col)
{
 double t, c, s, theta;
 int n = S.size1();
 //matrix<double> S(n,n);
 //S = S0;

 if (row == col) return;

 theta=(S(row,row)-S(col,col))/(2*S(row,col));
 if (theta < EPS)
 {
  t = 1/(abs(theta)+sqrt(theta*theta+1));
 }
 else {
  t = 1/abs(2*theta);
 }

 if (theta<0) t = -t;

 c = 1/sqrt(t*t + 1);
 s = c*t;
 matrix<double> R(n,n);
 R = identity_matrix<double>(n,n);
 R(row,row) = R(col,col) = c; R(row,col) = s; R(col,row) = -s;
 S = prod(S, trans(R));
 S = prod(R, S);
 U = prod(U, R);
}


int jacobi2(matrix<double> &S, vector<double> &e, matrix<double>  &U)
{
 int col, row;
 bool iterating = true;
 int n = S.size1();
 if (S.size2() != n)
 {
  return -1;
 }
 matrix<double> M(n,n);

 U = identity_matrix<double>(n,n);
 
 while(iterating)
 {
  M = S;
  abs(M);

  for (int k=0; k<n; k++)
  {
   M(k,k)=0;
  }

  findMax(M, row, col);
  if (row == col)
  {
   for (int i=0; i<n; i++) e(i) = S(i,i);
   return 0;
  }
  double Smax = S(row,col);
  zeroRowCol (S, U, row, col);
  if (Smax < EPS * norm_frobenius(S)) iterating = false;
 }

 for (int i=0; i<n; i++) e(i) = S(i,i);

 return 0;  
} 


int _tmain(int argc, _TCHAR* argv[])
{
 int n = 3;

 mmatrix<double> M(n,n);
 M = boost::numeric::ublas::identity_matrix<double>(3,3);
 M(0,1)=.5;
 M(1,0)=.5;

 std::cout << "The matrix is:\n" << M << std::endl;

 matrix<double> M(n,n);
 matrix<double> U(n,n);
 vector<double> e(n); 

 jacobi2(M, e, U);

 std::cout << M << std::endl;
 std::cout << e << std::endl;
 std::cout << U << std::endl;
}

Wednesday, June 27, 2012

The Fourier transform as a diagonalization


One of the benefits of using the Fourier transform of a function is that convolutions become multiplications. This is important when solving a differential equation with its Green's function. If the Green's function comes from a differential operator $D^*D$, where $D$ is a differential operator and $D^*$ is its adjunct, then the Green's function is not singular at the origin, and is continuous. It expands a function space called a reproducing kernel hilbert space, RKHS, and all functions in this space can be written as linear combinations of the Green's function evaluated on one argument, and the solution to the differential equation $D^*D u = y$ would be of that form. OK, don't digrees anymore... to the cheese...

In the Fourier domain we operate on frequencies $\omega$. For example, to attenuate the noise, we decrease the power in the high omegas, which accounts for a convolution (with a Gaussian, for example). If we see this linear operation as a matrix, the convolution operator that has one (let's say) dimensional Gaussians in its rows (in the time/space domain) becomes a diagonal in the Fourier domain.

The page popped up with much to follow on. In particular, I liked this paragraph
The moral of the story is that the Fourier Transform may be thought of as a change of basis.  The Fourier integral projects a function onto the basis functions of a new coordinate system whose basis functions are the complex exponentials.  In this new basis, the convolution operator is diagonal and everything is simple.  The convolution operator acts on each Fourier component independently by multiplying the component by an associated magnitude and phase.
In Matlab
C=[4 1 2 3; 3 4 1 2; 2 3 4 1; 1 2 3 4]

C =

     4     1     2     3
     3     4     1     2
     2     3     4     1
     1     2     3     4

 F=fft(C)

F =

  10.0000            10.0000            10.0000            10.0000        
   2.0000 - 2.0000i  -2.0000 - 2.0000i  -2.0000 + 2.0000i   2.0000 + 2.0000i
   2.0000            -2.0000             2.0000            -2.0000        
   2.0000 + 2.0000i  -2.0000 + 2.0000i  -2.0000 - 2.0000i   2.0000 - 2.0000i

  F*C*F'

ans =

  1.0e+003 *

   4.0000                  0                  0                  0          
        0             0.0640 - 0.0640i        0                  0          
        0                  0             0.0320                  0          
        0                  0                  0             0.0640 + 0.0640i