iconEuler Reference

Linear Algebra

Linear algebra and regression routines.

This file contains basic linear algebra routines, like solutions of linear systems, linear fit, polynomial fit, eigenvalues and eigenspaces, decompositions, singular values, norms, and special matrices.

It does also contain methods to solve systems by hand (pivotize).

Linear Systems

function id (n:index, m:index=none)

  Creates a nxn identity matrix.
  
  See: 
diag (Euler Core),
diag (Maxima Documentation),
eye (Linear Algebra)
function eye (n:positive integer,m:index=none)

  Creates a nxm matrix with 1 on the diagonal
  
  If m=none, then it creates an nxn matrix. This is a Matlab
  function.
  
  See: 
id (Linear Algebra)
function inv (A:complex)

  Computes the inverse of a matrix, using the Gauss algorithm.
  
  See: 
xinv (Exact Computation),
svdsolve (Linear Algebra),
xlgs (Exact Computation),
fit (Linear Algebra),
pinv (Linear Algebra)
function invert (A)

  returns A^{-1}
  Alias to inv.
  
  See: 
pinv (Linear Algebra)
function transpose (A)

  returns A'
  Alias for Maxima
function image (A:complex)

  Computes the image of the quadratic matrix A.
  
  Returns a matrix, containing a basis of the image of A in the
  columns.
  
  The image of a matrix is the space spanned by its columns. This
  functions uses an LU decomposition with the Gauss algorithm to
  determine a base of the columns. The function "svdimage" returns an
  orthogonal basis of the image with more effort.
  
  >A=redim(1:9,3,3)
  1                   2                   3
  4                   5                   6
  7                   8                   9
  >rank(A)
  2
  >image(A)
  0.0592348877759      0.118469775552
  0.236939551104       0.29617443888
  0.414644214431      0.473879102207
  >svdimage(A)
  -0.214837238368     -0.887230688346
  -0.520587389465     -0.249643952988
  -0.826337540561       0.38794278237
  >kernel(A)
  1
  -2
  1
  >svdkernel(A)
  0.408248290464
  -0.816496580928
  0.408248290464
  
  See: 
svdimage (Linear Algebra)
function kernel (A:complex)

  Computes the kernel of the matrix A.
  
  Returns a matrix M, containing a basis of the kernel of A. I.e.,
  A.M=0.
  
  This function uses an LU decomposition with the Gauss algorithm. The
  function "svdkernel" does the same with more effort, but it returns
  an orthogonal base of the kernel.
  
  Add eps=..., if A is almost regular. Use a larger epsilon in this
  case than the default epsilon().
  
  See: 
image (Linear Algebra),
image (Maxima Documentation),
svdkernel (Linear Algebra)
function rank (A:complex)

  Computes the rank of the matrix A.
  
  The rank is the dimension of the image, the space spanned by the
  columns of A.
  
  See: 
image (Linear Algebra),
image (Maxima Documentation)
function totalsum (A)

  Computes the sum of all elements of A.
function seidel (A:real, b:real column, x=0, om=1.2)

  Solve Ax=b using the Gauß-Seidel method for large, sparse A.
  
  A must be diagonal dominant, at least not have 0 on the diagonal.
  The method will converge for all positiv definit matrices. However,
  larger dominance in the diagonal speeds the convergence. For thin
  matrices use seidelX.
  
  om is a relaxation parameter between 1 and 2. The default is
  1.2, and works in many cases. Specify the accuracy with an optional
  parameter eps. The optional column vector x is start value.
  
  See: 
seidelX (Large and Sparse Matrices)
function cg (A:real, b:real column, x0:column=none, f:index=10)

  Conjugate Gradent method to solve Ax=b.
  
  This function is the function of choice for large matrices. There
  is a variant "cgX" for large, sparse matrices.
  
  Stops as soon as the error gets larger.
  
  See: 
cgX (Large and Sparse Matrices),
cpxfit (Large and Sparse Matrices)
function solvespecial (A:complex, b:complex)

  Special solution of Ax=b using the LU decomposition.
  
  Returns a vector, such that A.x=b, even for singular, or non-square
  matrices A.
  
  This function uses the Gauss algorithm to determine the LU
  decomposition. It then extracts a base from the columns, and
  uses this base to get a special solution. Works only, if the matrix
  has maximal rank, i.e., the rows are linear independent.
  
  For matrices, which have a solution, the functions "svdsolve" or
  "fit" provide the same functionality with more effort.
  
  See: 
svdsolve (Linear Algebra),
fit (Linear Algebra),
kernel (Linear Algebra)
function det (A:complex)

  Determinant of A.
  
  Uses the Gauss algorithm to compute the determinant of A.
function cholesky (A)

  Decompose a real matrix A, such that A=L.L'. Returns L.
  
  A must be a positive definite symmetric and at least 2x2 matrix.
  A is not checked for symmetry,
  
  See: 
LU (Linear Algebra),
lu (Euler Core)
function lsolve (L,b)

  Solve L.L'.x=b
  
  See: 
cholesky (Linear Algebra),
cholesky (Maxima Documentation)
function choleskyeigen (A)

  Iterates the Cholesky-iteration, until the diagonal converges.
  
  A must be positive definite symmetric and at least 2x2.
  
  See: 
cholesky (Maxima Documentation)
function pivotize (A:numerical, i:index, j:index)

  Make A[i,j]=1 and all other elements in column j equal to 0.
  
  The function is using elementary operations on the rows of A. Use
  this functions for the demonstration of the Gauss algorithm.
  
  The function modifies A, but returns A, so that it can be printed
  easily.
  
  See: 
normalizeRow (Linear Algebra),
swapRows (Linear Algebra),
gjstep (Linear Algebra)
function pivotizeRows (A:numerical, i:index, j:index)

  Alias to pivotize
  
function normalizeRow (A:numerical, i:index, j:index)

  Divide row j by A[i,j].
  
  The function is using elementary operations. Use this functions for
  the demontration of the Gauss algorithm.
  
  See: 
pivotizeRows (Linear Algebra),
swapRows (Linear Algebra)
function swapRows (A:numerical, i1:index, i2:index)

  Swap rows i1 and i2 of A.
  
  The function is using elementary operations. Use this functions for
  the demontration of the Gauss algorithm.
  
  See: 
pivotizeRows (Linear Algebra),
normalizeRow (Linear Algebra)

QR Decomposition

function qrdecomp (A:real)

  Decompose QA=R using Givens rotations.
  
  The function returns {R,Q,c}. Q is an orthogonal matrix, and R is
  an upper triangle matrix. c is a vector with ones in the basis
  columns of A.
  
  This function uses orthogonal transformations to compute Q and R.
  It works for non-square matrices A.
  
  see: givensqr, givensrot

Fit

function fitnormal (A:complex, b:complex)

  Computes x such that ||A.x-b||_2 is minimal using the normal equation.
  
  A is an nxm matrix, and b is a nx1 vector.
  
  The function uses the normal equation A'.A.x=A'.b. This works only,
  if A has full rank. E.g., if n>m, then the m columns of A must be
  linear independent. The function "fit" provides a faster and better
  solution. "svdsolve" has the additional advantage that it produces
  the solution with minimal norm. For sparse matrices use "cpxfit"
  
  For badly conditioned A, you should use svdsolve instead.
  
  A : real or complex matrix (nxm)
  b : column vector or matrix (mx1 or mxk)
  
  See: 
svdsolve (Linear Algebra),
cpxfit (Large and Sparse Matrices)
function fit (A:real, b:real)

  Computes x such that ||A.x-b||_2 is minimal using Givens rotations.
  
  A is a nxm matrix, and b is a nxk vector (usually k=1).
  
  This function works even if A has no full rank. Since it uses
  orthogonal transformations, it us also quite stable. The function
  "svdsolve" does the same, but minimizes the norm of the solution
  too.
  
  A : real matrix (nxm)
  b : real column vector or matrix  (mx1 or mxk)
  
  See: 
fitnormal (Linear Algebra),
svdsolve (Linear Algebra)
function divideinto (A: real, B: real)

  Divide A into B yielding fit(B',A')'
  
  This function is used in Matlab mode for A/B.
  
  See: 
fit (Linear Algebra)
function polyfit (x:complex vector, y:complex vector, ..
    n:nonnegative integer scalar, w:real vector=none)

  Fits a polynomial in L_2-norm to data.
  
  Given data x[i] and y[i], this function computes a polynomial p of
  given degree such that the sum over the squared
  errors (p(x[i])-y[1])^2 is minimimal.
  
  The fit uses the normal equation solved with an orthogonal
  decomposition.
  
  w : Optional weights (squared) for the
  

Norm

function norm (A:numerical, p:nonnegative real=2)

  Euclidean norm of the vector A.
  
  For p=2, the function computes the square root of the sum of
  squares of all elements of A, if A is a matrix. For p>0 it computes
  the p-norm. For p=0 it computes the sup-norm.
  
  Note: The function does not compute the matrix norm.
  
  See: 
maxrowsum (Linear Algebra)
function maxrowsum (A:numerical)

  Compute the maximal row sum of A
function crossproduct (v, w)

  Cross product between two 3x1 or 1x3 vectors.
  
  The function will work for two column or row vectors. The result is
  of the same form as v.

Eigenvalues

function eigenvalues (A:complex)

  Eigenvalues of A.
  
  Returns a complex vector of eigenvalues of A.
  
  This function computes the characteristic polynomial using
  orthogonal transformations to an upper diagonal matrix, with
  additional elements below the diagonal. It then computes the zeros
  of this polynomial with "polysolve".
  
  See: 
charpoly (Euler Core),
charpoly (Maxima Documentation),
polysolve (Euler Core),
svd (Euler Core),
jacobi (Euler Core),
jacobi (Linear Algebra),
jacobi (Maxima Documentation),
svdeigenvalues (Linear Algebra),
mxmeigenvalues (Maxima Functions for Euler)
function eigenspace (A:complex, l:complex vector)

  Returns the eigenspace of A to the eigenvaue l.
  
  The eigenspace is computed with kernel(A-l*id(n)). Since the
  eigenvalue may be inaccurate, the function tries several epsilons
  to find a non-zero kernel. A more stable function is
  "svdeigenspace".
  
  See: 
eigenvalues (Linear Algebra),
eigenvalues (Maxima Documentation),
svdeigenspace (Linear Algebra)
function eigen (A:complex)

  Eigenvectors and a basis of eigenvectors of A.
  
  Returns {l,x,ll}, where l is a vector of eigenvalues, x is a basis
  of eigenvectors, and ll is a vector of distinct eigenvalues.
  
  See: 
eigenvalues (Linear Algebra),
eigenvalues (Maxima Documentation),
eigenspace (Linear Algebra),
mxmeigenvectors (Maxima Functions for Euler)
function cjacobi (A:complex, eps=sqrt(epsilon))

  Jacobi method to compute the eigenvalus of a Hermitian matrix A.
  
  Returns a vector of eigenvalues.
  
  This function computes the eigenvalues of a real or complex matrix
  A, using an iteration with orthogonal transformations. It works
  only for Hermitian matrices A, i.e., A=conj(A').
  
  See: 
jacobi (Euler Core),
jacobi (Maxima Documentation)
function jacobi (A:real)

  Jacobi method to compute the eigenvalus of a symmetric matrix A.
  
  Returns a vector of eigenvalues.
  
  This function computes the eigenvalues of a real matrix A, using an
  iteration with orthogonal transformations. It works only for
  symmetric matrices A, i.e., A=A'.
  
  See: 
cjacobi (Linear Algebra)

Hilbert Matrix

hilbertfactor:=3*3*3*5*5*7*11*13*17*19*23*29;
function hilbert (n:index, f=hilbertfactor)

  Scaled nxn Hilbert matrix.
  
  The Hilbert matrix is the matrix (1/(i+j-1)). This function
  multplies all entries with a factor, so thatthe matrix can
  accurately be stored up to the 30x30 Hilbert matrix.
function hilb (n:integer, f=hilbertfactor)

  Scaled nxn Hilbert matrix.
  
  Alias for hilbert(n).
  
  See: 
hilbert (Linear Algebra)

Matrix Functions

function power (A:numerical, n:integer scalar)

  A^n for integer n for square matrices A.
  
  The function uses a binary recursive algorithm. For n<0, it will
  compute the inverse matrix first. For n<0, the matrix must be real
  or complex, not interval.
  
  For scalar values, use the faster x^n. Moreover, power(x,n) will
  not map to the elements of x or n.
  
  See: 
matrixpower (Linear Algebra)
function matrixpospower (A:numerical, n:natural)

  Returns A^n for natural n>=0 and a matrix A.
  
  Alias for power(A,n)
function matrixpower (A:complex, n:integer scalar)

  Returns A^n for integer n and a matrix A.
  
  Alias for power(A,n)
function matrixfunction (f:string, A:complex)

  Evaluates the function f for a matrix A.
  
  The function f must be an analytic function in the spectrum of A.
  A must have a basis of eigenvectors. This function will decompose A
  into M.D.inv(M), where D is a diagonal matrix, and apply f to the
  diagonal of D.
  
  See: 
matrixpower (Linear Algebra)
function LU (A)

  Returns {L,R,P} such that L.R=P.A
  
  This function uses the Gauß-Algorithm in lr(A) to find the LR
  decomposition of a regular square matrix A. L is a lower triangle
  matrix with 1 on the diagonal, R an upper right triangle matrix,
  and P a permutation matrix.
  
  See: 
lu (Euler Core)
function echelon (A)

  Forms A into echelon form using the Gauß algorithm
  
  See: 
LU (Linear Algebra),
lu (Euler Core)

Filter

function filter (b:complex vector, a:complex vector, ..
    x:complex vector, y:complex vector=none, zeros:integer=true)

  Apply a filter with a and b to start values x and y
  
  Computes y[n] recursively using
  
  a[1]*y[n]+...+a[m]*y[n-m+1]=b[1]*x[n]+...+b[k]*x[n-k+1]
  
  a[1] must not be zero.
  
  The start values on the right hand side are either 0 or the first k
  values of x, depending on the flag zeros (true by default). The
  start values on the left hand side are either 0 (the default), or
  the values of the parameter y.
  
  The size of a must exceed the size of x by 1.
  All input vectors must be real or complex row vectors or scalars.
  
  See: 
fold (Euler Core),
fftfold (Numerical Algorithms)

Singular Values

The singular value decomposition of Euler is based on the builtin function svd. It is used to compute an orthogonal basis of the kernel and the image of a matrix, the condition of a matrix, or the pseudo-inverse.

function svdkernel (A:real)

  Computes the kernel of the quadratic matrix A
  
  This function is using the singular value decomposition, and works
  for real matrices only.
  
  Returns an orthogonal basis of the kernel as columns of a matrix.
  
function svdimage (A:real)

  Computes the image of the quadratic matrix A
  
  This function is using the singular value decomposition, and works
  for real matrices only.
  
  Returns an orthogonal basis of the image as columns of a matrix.
  
  See: 
kernel (Linear Algebra),
image (Maxima Documentation),
svdkernel (Linear Algebra)
function svdcondition (A:real)

  Condition number based on a singular value decompostion of A
  
  Returns the qutient of the largest singular value divided by the
  smallest. 0 means singularity.
  
  A : real matrix
function svddet (A:real)

  Determinant based on a singular value decomposition of A
  
  A : real matrix
  
function svdeigenvalues (A:real)

  Eigenvalues or singular values of A
  
  For a symmetric A, this returns the eigenvalues of A For a
  non-symmetric A, the singular values.
  
  A : real matrix
  
  See: 
eigenvalues (Maxima Documentation)
function svdeigenspace (A:real,l:real)

  Returns the eigenspace of A to the eigenvaue l
  
function svdsolve (A:real,b:real)

  Minimize |A.x-b| with smallest norm for x
  
  The singular value decomposition is one way to solve the fit
  problem. It has the advantage, that the result will be the result
  with smallest norm, if there is more than one solution.
  
  A : real matrix
  
  See: 
fit (Linear Algebra)
function svdinv (A:real)

  The pseudo-inverse of A.
  
  The pseudo-inverse B of a matrix solves the fit problem to minimize
  |Ax-b| by x=B.b. It is computed in this function using an svd
  decomposition.
  
function pinv (A:real)

  The pseudo-inverse of A.
  
  Alias to svdinv
  
  See: 
svdinv (Linear Algebra)
function ginv (A:real)

  Gilbert inverse of a matrix A
  
  This inverse has the property A.G.A=A

Gauß-Jordan Scheme

Some functions for the demonstration of the Gauß-Jordan algorithm. These functions can be used for linear systems or for the simplex algorithm.

defaultgjdigits:=3;
defaultgjlength:=10;
function gjprint (A : real, n : positive integer vector,
    v : string vector, digits:nonnegative integer=none,
    length:index=none)

  Print the matrix A in Gauß-Jordan form.
  
  The scheme A is printed in the form
  
  x  y ...
  a  1  2 ...
  b  3  4 ...
  ...
  
  n : index vector with a permumatin of the variables
  v : variable names with the column variables first
  digits : number of digits for each number
  length : output length of each number
  
  File: gaussjordan
  
  See: 
gjstep (Linear Algebra)
function gjstep (A:numerical, i:index, j:index,
    n:positive integer vector=none,
    v:string vector=none, digits:nonnegative integer=none,
    length:index=none)

  Make A[i,j]=1 and all other elements in column j equal to 0.
  
  The function is using elementary operations on the rows of A. Use
  this functions for the demonstration of the Gauss algorithm.
  
  The function modifies A and n.
  
  n : A row vector of indices, the function assumes the Gauß-Jordan
  form. n contains the indices of the variables, first the indices of
  the variables in the rows, and then the indices of the variables in
  the columns. The row variable i is exchanged with the column
  variable j. This is the same as makeing the j-th column to an
  canonical vector e[i], and inserting in the j-th column the result
  of this operation applied to e[i].
  
  If n is present, you can add a vector v of strings, which contains
  the names of the variables. The problem will then be printed using
  gjprint.
  
  File: gaussjordan
  
  See: 
pivotize (Linear Algebra),
gjprint (Linear Algebra)
function gjsolution (M : real, n : positive integer vector)

  Read the values of the solution from the Gauß-Jordan scheme.
  
  For this, we assume that the last column contains the values, and
  the variables in the top row are set to zero.
  
  See: 
gjstep (Linear Algebra)
function gjvars (M : real, simplex=false)

  Generate variable names for M
  
  simplex : If true the last row variable will be named ""
  and the last column will name will be missing.
  
  Returns {v,n}
  where n=1:(rows+cols), v=["s1"...,"x1",...]
  
  See: 
gjprint (Linear Algebra)
function gjaddcondition (M : real,
    n:positive integer vector, v:string vector,
    line:real vector, varname:string)

  Add a condition to the Gauß-Jordan Form
  
  The line for the condition (left side <= right side) and the
  variable name must be provided. The condition uses the current top
  row variables.
  
  Return {M,n,v}

Documentation Homepage