Tuesday, October 2, 2012

Machine Learning video lectures (Beginners)

The first one is part of the PyData workshop held in March this year. This lecture is oriented to beginners, so if you already know the (very) basics, then it won't be of interest to you. It is nice, however, since you can see all the advanced capabilities of Scikits Learn while watching something you already know, such as linear, polynomial and Gaussian SVMs, logistic regression and naive bayes.


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;
}