Monday, November 14, 2011

Scilab - Submatrix Operations

Getting up to speed with sub-matrix operations in Scilab

For a beginner, sub-matrix operations appear a little baffling. You may need time to master the techniques of sub-matrix operations, but once you do, they are nothing short of amazing.

Let us first see why we need them. Sub-matrix operations can do two things:
  1. Extract parts of an existing matrix, and can stored in a new matrix for subsequent use or be used immediately in an expression.
  2. Replace parts of an existing matrix with another matrix. Of course, size of the sub-matix that you are replacing must be the same size as the matrix you want to place in it.
As with most cases, it is best explained with an example. Let us say we have a matrix of size 5x4 as shown below:
-->x = [1 2 3 4 5; 6 7 8 9 10; 11 12 13 14 15; 16 17 18 19 20]
x  =
   1.    2.    3.    4.    5.
   6.    7.    8.    9.   10.
  11.   12.   13.   14.   15.
  16.   17.   18.   19.   20.
Let us now extract elements of 'x' from rows 2 to 4 and columns 3 to 5, and place the sub-matrix of 'x' so extracted in a new matrix 'y'. This accomplished as follows:
-->y = x(2:4, 3:5)
y  =
   8.    9.   10.
  13.   14.   15.
  18.   19.   20.
Let us demystify this mumbo-jumbo:
  1. Within parentheses, the numbers 2:4 and 3:5 are called a range, another way of saying numbers starting from 2, up to and including 4, at an increment of 1 (or from 3 up to and including 5, at an increment of 1).
  2. 2:4 indicate the range of rows of 'x' that are to be extracted. 3:5 indicate the range of columns of 'x' to be extracted. Thus, the sub-matrix being extracted is of size 3x3 (rows 2, 3 and 4 and columns 3, 4 and 5).
  3. This portion of 'x' is copied (not literally extracted, but only copied) and placed in a matrix 'y'. If there are existing elements in 'y', they are discarded. If 'y' does not exist, it is created.
Actually, it is a straight forward extension of the concept of range. Thus, the following are possible:
  1. Row and/or column ranges can have an increment other than 1, as long as it is an integer. For example using an increment of 2 extracts alternate rows and/or columns.
  2. Row and/or columns ranges can be decreasing instead of increasing, and decrement can be 1 or any other negative value, as long as it is an integer. Range 5:-1:3 is 5, 4, 3. Range 4:-1:2 is 4, 3, 2.
When ranges are used in sub-matrix extraction, the following additional riles apply:
  1. The number of the last row/column can be represented by the dollar symbol $.
  2. While representing the range of rows/columns, if both the start and end value are left out, they are assumed to be 1 and $, respectively.
  3. It acceptable to specify both the start and end values of a range or leave out both. Leaving out only one of them results in an error.
So here are a few tricks you can decipher on your own:
  1. y = x(1:2:$, 1:2:$) copies odd rows and odd columns of 'x' into 'y'.
  2. y = x(2:2:$, 2:2:$) copies even rows and even columns of 'x' into 'y'.
  3. y = x(:, :) is the same as y = x
  4. y = x($:-1:1, :) copies rows of x in reverse order while not altering the columns. The last row of 'x' will be the first row of 'y', last but one row of 'x' will be the second row of 'y', and so on. The rows of 'x' are inverted and placed in 'y'.
  5. y = x($:-1:1, $:-1:1) reverses both rows and columns.
  6. y = x(1:$-1, 1:$-1) copies rows from 1 to the last but one and columns starting from 1 to last but one. Notice, $ is used just as you would use a variable in a programming language.
  7. y = x(1:, :) is an error because in specifying the range of rows, the end value is not mentioned. Column range specification is acceptable, and actually represents all columns.
Whatever you do with copying/extraction works equally well the other way around.
-->x(2:4, 3:5) = zeros(3, 3)
Here, on the right hand side of the expression we create a matrix of size 3x3 containing all zeros, and place this matrix into the sub-matrix of 'x' between rows 2 to 4 and columns 3 to 5. The matrices on either side of the assignment operator must be of the same size. If you don't want to bother calculating the size of the zero matrix, you can leave it to Scilab to do all the hard work. Thus, the following commands work:
-->x(2:4, 3:5) = zeros()
-->x(1:4, 1:2) = ones()

Well, that is neat thinking by people who designed Scilab programming language. But it is not unique. Matlab(R) and GNU Octave do almost the same thing. SO you learn sub-matrix operations in Scilab, you have a head start if you ever have to use either of them.

2 comments:

  1. also note that M(:) transform any matrix into a column vector :

    -->M=rand(3,3)
    M =

    0.2113249 0.3303271 0.8497452
    0.7560439 0.6653811 0.6857310
    0.0002211 0.6283918 0.8782165

    -->M(:)
    ans =

    0.2113249
    0.7560439
    0.0002211
    0.3303271
    0.6653811
    0.6283918
    0.8497452
    0.6857310
    0.8782165

    this is usefuill to study linear application over a matrix space (as M --> A*M+M*A)

    ReplyDelete
  2. Philippe, thanks for pointing it out. In fact, Scilab allows you to treat a matrix (2D) as a vector (1D), as long as you remember that elements are indexed column-wise. That is, matrix indexes are mapped to vector indexes in the sequence: (1,1) is 1, (2,1) is 2, ..., (m,1) is m+1. In general, for a matrix with 'm' rows and 'n' columns, the element (i, j) of the matrix is element (j-1)*m+i of the vector.

    In fact, this is the internal storage scheme used in Fortran (assuming start index as 1). C uses the row-wise scheme (and with start index as zero, the mapping function is (i*n+j))

    ReplyDelete

Your comments are ciritical input for the sustainance and improvement of this blog.