Skip to content
jcanny edited this page Jan 29, 2013 · 3 revisions

Table of Contents

Element-wise Functions

Currently, the only functions in this category are max and min (most element-wise functions are listed under scientific functions later, and element-wise operators are list under Matrix algebra). max(A,B) for equal-sized A and B computes the element-wise max of the matrices. Just as for matrix operators, matrix functions also work in these cases:

  • A or B is a scalar, in which case that value is applied against all values in the other matrix.
  • A or B is a 1 x 1 matrix, in which case the value is treated as a scalar as above.
  • A or B is a vector whose length matches the corresponding dimension of the other matrix. In that case, each value in the vector is applied against all the elements in the other matrix in the corresponding row or column.

Reducing Functions

Reducing functions apply a binary operator cumulatively to rows or columns of the input matrix, and produce either a vector output (standard reducing functions) or a same-size matrix (cumulative reducing functions). The standard reducing functions are

  • maxi(A, n), mini(A, n), sum(A, n), mean(A, n), variance(A, n), sdev(A, n)
Each of these reduces along the dimension specified by n. e.g. maxi(A,1) computes a row vector which is the max of each column, while maxi(A,2) computes a column vector whose values are the row-wise maximum values of A.

All of these functions can be called without a dimension argument, which is a convenience form intended for vector arguments. So e.g. sum(A) for a vector A reduces along the long dimension of A. Calling sum(A) without the dimension argument for a non-vector A will reduce along dimension 1.

The cumulative reducing functions right now are:

  • cumsum(A, n)
which computes the cumulative sum along dimension n, and returns every partial sum. cumsum(A) for vector A works along the long dimension of A.

The following function constructs dense matrices from tuples of indices.

  • accum(inds:IMat, vals:DMat, nrows, ncols)
sums all the values for a given index value, and stores the result in that location. The inds matrix can be either a vector of single indices, or a 2 x n matrix whose columns are row and column indices respectively. If called without dimension arguments: accum(inds, vals) the matrix is sized to fit the largest row and column indices.

Sorting Functions

The sort functions sorts columns of a matrix in ascending order. sort2(A) returns both the a sorted version of A and a permutation matrix P such that for each i, B(?,i) = A(P(?,i),i). Or if we add column offsets to P, POFF, then we have that B=A(P+POFF).

scala> a
res16: BIDMat.DMat =
   0.14771   0.53622   0.85483  0.042617
   0.96534   0.41199   0.63428  0.014676
   0.66535  0.029090   0.85247   0.61777
   0.38885   0.50339   0.81829   0.77564

scala> val (b, p) = sort2(a)
b: BIDMat.DMat =
   0.14771  0.029090   0.63428  0.014676
   0.38885   0.41199   0.81829  0.042617
   0.66535   0.50339   0.85247   0.61777
   0.96534   0.53622   0.85483   0.77564

p: BIDMat.IMat =
   0   2   1   1
   3   1   3   0
   2   3   2   2
   1   0   0   3

scala> val poff = iones(4,1)*irow(0 to 12 by 4)
poff: BIDMat.IMat =
   0   4   8  12
   0   4   8  12
   0   4   8  12
   0   4   8  12

scala> a(p + poff)
res21: BIDMat.DMat =
   0.14771  0.029090   0.63428  0.014676
   0.38885   0.41199   0.81829  0.042617
   0.66535   0.50339   0.85247   0.61777
   0.96534   0.53622   0.85483   0.77564

scala> b
res22: BIDMat.DMat =
   0.14771  0.029090   0.63428  0.014676
   0.38885   0.41199   0.81829  0.042617
   0.66535   0.50339   0.85247   0.61777
   0.96534   0.53622   0.85483   0.77564

sort and sort2 always sort in ascending order. For descending sorts, there are functions sortdown and sortdown2.

To sort each row, use sort(A,2) or sortdown(A,2).

To sort rows of a matrix (in lexicographic order), there are two functions sortlex(A) and sortrows(A). sortlex returns only the permutation matrix for the sort. sortrows returns both the permutation matrix and the matrix itself.

scala> a
a: BIDMat.DMat =
   0.083308  0.0053793    0.60398
    0.75904    0.91294    0.43385
    0.53446    0.91638    0.59250
    0.39092    0.65965   0.020080

scala> val (aa, ii) = sortrows(a)
aa: BIDMat.DMat =
   0.083308  0.0053793    0.60398
    0.39092    0.65965   0.020080
    0.53446    0.91638    0.59250
    0.75904    0.91294    0.43385

ii: BIDMat.IMat =
   0
   3
   2
   1

scala> a(ii,?)
res26: BIDMat.DMat =
   0.083308  0.0053793    0.60398
    0.39092    0.65965   0.020080
    0.53446    0.91638    0.59250
    0.75904    0.91294    0.43385

Finally, for descending sorts there are functions sortlexdown and sortrowsdown.

Set Functions

unique and uniquerows

Matrix Inverse

The function inv(A) returns the inverse of the square matrix A. Matrix inversion is an O(n^3) operation.

Eigenvalues

There are two eigenvalue functions for symmetric matrices. Both can return two values: a column vector of eigenvalues, and a square matrix whose columns are the corresponding eigenvectors. The eig function takes two arguments, the matrix A, and a boolean flag which determines whether the eigenvectors are returned. eig uses reduction to tridiagonal form followed by calculation of eigenvalues. The second function is feig which uses a faster, divide-and-conquer algorithm for eigenvalues and eigenvectors.

scala> aa
aa: BIDMat.DMat =
  2.9260  1.6648  2.5727  2.2300
  1.6648  1.3799  1.6853  1.4514
  2.5727  1.6853  3.1019  2.3740
  2.2300  1.4514  2.3740  2.2825

scala> val (v, m)= eig(aa, true)
v: BIDMat.DMat =
  0.27301
  0.31132
  0.45670
   8.6492

m: BIDMat.DMat =
  -0.13731   0.40833  -0.71388   0.55207
  -0.12716  -0.90024  -0.21440   0.35697
  -0.50598   0.14970   0.62654   0.57360
   0.84200  0.020595   0.22771   0.48863

scala> aa * m /@ m
res31: BIDMat.DMat =
  0.27301  0.31132  0.45670   8.6492
  0.27301  0.31132  0.45670   8.6492
  0.27301  0.31132  0.45670   8.6492
  0.27301  0.31132  0.45670   8.6492

Both algorithms take O(n^3) steps, and feig is about twice as fast as eig when retrieving the eignvectors. feig should be used only with positive definite matrices. eig requires about twice as many flops (3 n^3 vs 4/3 n^3) as matrix inversion, but the calculation is less uniform and the MFlop/s achieved is lower.

Cholesky Decomposition

The function chol computes a Cholesky factorization A=L*L^T of a symmetric matrix. The factor is returned a a lower-triangular matrix.

scala> aa
res33: BIDMat.DMat =
  2.9260  1.6648  2.5727  2.2300
  1.6648  1.3799  1.6853  1.4514
  2.5727  1.6853  3.1019  2.3740
  2.2300  1.4514  2.3740  2.2825

scala> val b = chol(aa)
b: BIDMat.DMat =
   1.7105        0        0        0
  0.97325  0.65778        0        0
   1.5040  0.33675  0.85233        0
   1.3037  0.27752  0.37518  0.60420

scala> b * b.t
res34: BIDMat.DMat =
  2.9260  1.6648  2.5727  2.2300
  1.6648  1.3799  1.6853  1.4514
  2.5727  1.6853  3.1019  2.3740
  2.2300  1.4514  2.3740  2.2825

Cholesky factorization is an O(n^3) calculation, and is about twice as fast as matrix inversion.

Triangular Solve

The trisolve function solves triangular systems of the form T x = y for x. trisolve accepts a string argument that specifies the type of calculation. The mode string has three letters:

  • Char1 is "U" or "L" for upper or lower triangular matrices.
  • Char2 = "N", "T" or "C" for A not-transposed, transposed or conjugate respectively.
  • Char3 = "N" or "U" whether the leading diagonal is non-unit "N" or unit "U" respectively.
scala> b
res35: BIDMat.DMat =
   1.7105        0        0        0
  0.97325  0.65778        0        0
   1.5040  0.33675  0.85233        0
   1.3037  0.27752  0.37518  0.60420

scala> val y = rand(4,1)
y: BIDMat.DMat =
  0.047913
   0.65338
   0.58416
   0.26086

scala> val x = trisolve(b, y, "LNN")
x: BIDMat.DMat =
  0.028011
   0.95187
   0.25987
  -0.22728

scala> b*x
res36: BIDMat.DMat =
  0.047913
   0.65338
   0.58416
   0.26086

Triangular matrix solving is an O(n^2) operation for vectors x and y, or O(n^2 k) for k x n matrices x, y.