Skip to content

Linear Solvers and Eigenvalues

Jump to bottom
jcanny edited this page May 30, 2014 · 12 revisions

Table of Contents

Matrix Inverse

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

Eigenvalues

Real Symmetric Matrices

There are two eigenvalue functions for real symmetric matrices seig and feig. Both can return two values: a column vector of eigenvalues, and a square matrix whose columns are the corresponding eigenvectors. The seig function takes two arguments, the matrix A, and a boolean flag which determines whether the eigenvectors are returned. seig 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. Both functions work on both single precision (FMat) and double precision (DMat) arguments. 

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 seig when retrieving the eignvectors. feig should be used only with positive definite matrices. seig 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.

General Real Matrices

There is a third eigenvalue routine which can deal with non-symmetric matrices. This routine is called geig and returns complex results.

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.

Clone this wiki locally