Because best way to understand is explaining.

Monday, June 4, 2012

Gentle Intro to Octave or Matlab

I began using Octave for homework assignments from the online Machine Learning class. Having worked with languages like Python, Groovy and JavaScript I never expected a system designed for numerical computations to include such a complete and unique programming language. But it does and I can't resist sharing some examples.

There are two important things one should know about Octave (or Matlab as Octave is usually portable to Matlab):
  • Octave is a high level language just like Python or Groovy
  • Using Octave without matrices or vectors is like using Java without objects
Just these by themselves are worth a whole book on Octave but instead I go on with few cool  examples (leaving the book for later :-).

Matrices and Vectors

Creating vector or matrix in Octave is simple:
octave:> A = [1 2 3; 4 5 6; 7 8 9]
A =

   1   2   3
   4   5   6
   7   8   9
defines 3x3 matrix of integers.

Use special functions to define special matrices, e.g. identity:
octave:> I = eye(3)
I =

Diagonal Matrix

   1   0   0
   0   1   0
   0   0   1

all zeros:

octave:> allZeros = zeros(2,4)
allZeros =

   0   0   0   0
   0   0   0   0
vector (number of columns is 1) of all ones:

octave:> allOnes = ones(3,1)
allOnes =

   1
   1
   1
or matrix with random values:

octave:> X = rand(3, 5)
X =

   0.400801   0.091597   0.951333   0.063074   0.018309
   0.690633   0.194094   0.417911   0.658953   0.624323
   0.848887   0.696741   0.213559   0.363656   0.632738
And finally getting a vector of values from 1 to N (row vector):
octave:> 1:N
ans =

    1    2    3    4    5    6    7    8    9   10
 and column (vector above transposed):
octave:> (1:N)'
ans =

    1
    2
    3
    4
    5
    6
    7
    8
    9
   10

To stir things up a bit I make the following claim:
In Octave for any given problem there is higher than 50% chance that using matrices alone solves the problem with less code and more efficiently than when using loop and condition statements. 
Being a high level language Octave has control statements if, switch, loops for and while but using them in Octave is often your second choice. The reason are many matrix operators and functions Octave offers may accomplish a task without ever invoking a single control statement in a fraction of time.

Suppose you have a matrix X and you need to insert a column of 1s in front. Then I just concatenate a vector of 1s of proper size and X:

 X = [ ones(size(X, 1), 1),  X ];

What just happened: function size returned 1st dimension of array X (number of rows); then function ones generated a vector (2d dimension is 1) of 1s, and finally we concatenated column and X. But this is only beginning.

Matrix magic

This example illustrates why and how things may work out better without control statements in Octave. Suppose I have a row vector (we may call it also an array but ultimately it is a single row matrix) of numbers from 1 to 10:
octave:> y = randperm(10)
y =

    4    3    1    8    6   10    9    2    7    5

octave:> y = repmat(y, 1, 10);

Function randperm returned a row vector containing random permutation of numbers 1 through 10. Then I used function repmat to make vector y 10 times its original length by repeating it 10 times (which results in 1 by 100 matrix y).

Now, to the problem we will solve. Let's say each number from 1 to 10 corresponds to 10-dimensional vector of 0s and 1 with all of its elements set to 0 except the position of the number. Thus, 1 corresponds to vector (1 0 0 0 0 0 0 0 0 0), 2 corresponds to vector (0 1 0 0 0 0 0 0 0 0) and so on until 10 that corresponds to (0 0 0 0 0 0 0 0 0 1). Then given some vector y like above (where each element is a number from 1 to 10) we want to produce series of vectors that correspond to elements in y.

If you never worked with Octave before then solution may amaze you, if you did then this might be your normal routine:

A = eye(10)

A =

Diagonal Matrix

   1   0   0   0   0   0   0   0   0   0
   0   1   0   0   0   0   0   0   0   0
   0   0   1   0   0   0   0   0   0   0
   0   0   0   1   0   0   0   0   0   0
   0   0   0   0   1   0   0   0   0   0
   0   0   0   0   0   1   0   0   0   0
   0   0   0   0   0   0   1   0   0   0
   0   0   0   0   0   0   0   1   0   0
   0   0   0   0   0   0   0   0   1   0
   0   0   0   0   0   0   0   0   0   1

result = A(:, y);

What just happened: first we created an identity matrix A of size 10 - note that it consists of exactly 10 vectors we are mapping to and each is in right column position. Now, we just plug our original vector y into column index of matrix A. This will extract elements from A: all for rows and precisely right columns. Thus A(:, 2) gives us matrix which is 2nd column of A, A(:, 2:4) gives us matrix with columns of A from 2 to 4, same is accomplished with A(:, [2:4]), and A(:, [1 4 9]) selects 1st, 4th and 9th columns of A. Finally, we can plug an arbitrary vector in column index - in our case vector y just what we need and result becomes 10 by 100 matrix where each column corresponds to element of y.

As unconventional as it may sound I keep thinking of this solution in terms of ranges when indexing arrays. Indexing array could be via an integer, range, or array of numbers. The latter is just a vector and that is all to it.

No comments: