iconEuler Reference

Linear Algebra

Content

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).

For an introduction and the core functions see the following links.

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 diagmatrix (v:vector, k:integer=0)
  Create a diagonal square matrix with v on the k-th diagonal
  
  The size of the square matrix is chosen, such that the vector v
  fits into the matrix.
  
  diagmatrix(1:3)
  diag(3,3,0,1:3) // the same
  diagmatrix(1:3,1)
  diag(4,4,1,1:3) // the same
  
  See: 
diag (Euler Core),
diag (Maxima Documentation)
function transpose (A)
  returns A'
  Alias for Maxima
function image (A:complex, eps=none)
  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, eps=none)
  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)

Matrix Functions

function overwrite max 
  max(A) or max(x,y,...)
  
  Returns the maximum of the rows of A for one argument,
  or the maximum of x,y,z,... for several arguments.
  
  See: 
max (Maxima Documentation),
extrema (Linear Algebra)
function overwrite min 
  min(A) or min(x,y,...)
  
  Returns the minimum of the rows of A for one argument,
  or the minimum of x,y,z,... for several arguments.
  
  See: 
min (Maxima Documentation),
extrema (Linear Algebra)
function comment extrema (A)
  Extremal values [min,imin,max,imax] in each row of A
  
  >x=random(5)
  [0.866587,  0.536137,  0.493453,  0.601344,  0.659461]
  >extrema(x)
  [0.493453,  3,  0.866587,  1]
  
  See: 
min (Linear Algebra),
min (Maxima Documentation),
max (Linear Algebra),
max (Maxima Documentation)
function comment totalmax (A)
  Maximal entry in the matrix A
function comment totalmin (A)
  Minimal entry in the matrix A
function comment flipx (A)
  flips a matrix vertically
function comment flipy (A)
  flips a matrix horizontally
function flipud (A:numerical)
  flips a matrix vertically.
function fliplr (A:numerical)
  flips a matrix horizontally.
function comment rotleft (A)
  rotates the rows of A to the left
function comment rotright (A)
  rotates the rows of A to the right
function comment shiftleft (A)
  shifts the rows of A to the left, last column to 0
function comment rotright (A)
  shifts the rows of A to the right, first column to 0
function rot90 (A:numerical)
  rotates a matrix counterclockwise
function rot (A:numerical, k:integer=1)
  rotates a matrix counterclockwise k times 90 degrees
function overwrite drop (A, i, cols=0)
  Drop rows i from A, or elements i from row vector
  
  cols: If true, the columns i are dropped from the matrix A.
  
  Overwrites the built-in function drop, which works for real row
  vectors only.

For sums and products of vectors and rows of matrices have a look at the core functions.

function totalsum (A)
  Computes the sum of all elements of A.
function flatten (A)
  Reshapes A to a vector

Solving Linear Systems

The simple Gauss algorithm is implemented in the A\b operator in Euler and in an algorithm from AlgLib.

The Gauss-Seidel method and the conjugate gradient methods are implemented for full and sparse matrices. Moreover, there are decompositions in LU or LL', and routines to solve systems based on the decompositions.

Euler has functions for sparse systems, exact functions, and interval routines for linear systems.

function overwrite alsolve (A:real, B:real,
    exact=0, check=1)
  Solve A.X=B
  
  A must be square, and B can have more than one column. The function
  uses an implementation of the Gauss algorithm from AlgLib. The
  result is of the same order as the Euler operator A\B.
  
  exacter : If set, refinements are tried to get a better solution.
  check : If set, low condition numbers issue an error exception.
  
  Returns {X,info,condition}
  info : If -3, the matrix was found to be singular
  condition : If close to 0, the matrix should be singular (inverse
  condition)
  
  >A=normal(200,200); b=sum(A);
  >longest totalmax(abs(alsolve(A,b)-1))
  1.321165399303936e-014
  >longest totalmax(abs(A\b-1))
  3.708144902248023e-014
  
function inv (A:complex)
  Computes the inverse of A.
  
  See: 
xinv (Numerical Algorithms),
svdsolve (Linear Algebra),
xlgs (Numerical Algorithms),
fit (Linear Algebra),
pinv (Linear Algebra)
function invert (A)
  Computes the inverse of A.
  
  Alias to inv.
function comment lu (A)
  LU-decomposition of A.
  
  The result is {B,r,c,det}.
  
  B is the result of the Gauss algorithm with the factors written
  below the diagonal, r is the index rearrangement of the rows of A,
  c is 1 for each linear independent column and det is the
  determinant of A.
  
  There is a utility function LU, which computes the LU-decomposiotn
  of a matrix (L,R,P such that LR=PA).
  
  To get a true LU-decomposition using lu for a non-singular square
  A, take the lower part of B (plus the identity-matrix) as L and the
  upper part as U. Then A[r] is L.U, i.e. for quadratic A, A[r] is
  
  >(band(B[r],-n,-1)+id(n)).band(B[r],0,n)
  
  To solve A.x=b for quadratic A and an x quickly, use
  lusolve(B[r],b[r]).
  
  >A=normal(20,20); b=sum(A);
  >{B,r,c,det}=lu(A); det,
  -2.61073e+008
  >longest totalmax(abs((band(B[r],-20,-1)+id(20)).band(B[r],0,20)-A[r]))
  2.220446049250313e-015
  >longest totalmax(abs(lusolve(B[r],b[r])-1))
  1.06581410364015e-014
  
  See: 
lusolve (Linear Algebra)
function comment lusolve (R,b)
  Solve linear system LUx=b.
  
  L has ones on the diagonal. It is stored in the lower left part of
  R. U is stored in the upper right part of R. This is the way that
  lu returns the matrices.
  
  >R=[1,2,3;0,4,5;0,0,6]
  1          2          3
  0          4          5
  0          0          6
  >L=[0;1;1,2,0]
  0          0          0
  1          0          0
  1          2          0
  >A=(L+id(3)).R
  1          2          3
  1          6          8
  1         10         19
  >A\sum(A)
  1
  1
  1
  >lusolve(R+L,sum(A))
  1
  1
  1
  
function seidel (A:real, b:real column, x=0, om=1.2, eps=none)
  Solve Ax=b using the Gauss-Seidel method.
  
  A must be diagonal dominant, at least not have 0 on the diagonal.
  The method will converge for all positive definite matrices.
  However, larger dominance in the diagonal speeds up the
  convergence. For thin matrices use seidelX().
  
  om :
  This is a relaxation parameter between 1 and 2. The default is 1.2,
  and works in many cases.
  
  x : start value
  
  You can specify the accuracy with an optional parameter eps.
  
  See: 
seidelX (Numerical Algorithms)
function cg (A:real, b:real column, x0:column=none,
    f:index=10, eps=none)
  Conjugate gradient 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 (Numerical Algorithms),
cpxfit (Numerical Algorithms)
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, eps=none)
  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 (Linear Algebra)
function lsolve (L,b)
  Solve L.L'.x=b
  
  See: 
cholesky (Linear Algebra),
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)
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)
function comment givensrot (j,i1,i2,A,B)
  One Givens rotation of A and B
  
  Returns (X=GA,Y=GB), such that X[l2,j]=0, and G is an orthogonal
  Matrix, changing only the rows i1 and i2 of A and B.
  
  >shortformat; A=random(3,3),
  0.77473     0.89357      0.5985
  0.70968     0.36746     0.80751
  0.59486     0.67635     0.26406
  >B=id(3)
  1           0           0
  0           1           0
  0           0           1
  >{X,Y}=givensrot(1,1,2,A,B); X
  -1.0506    -0.90712    -0.98678
  0     0.33263    -0.19118
  0.59486     0.67635     0.26406
  >Y.Y'
  1           0           0
  0           1           0
  0           0           1
  
  See: 
givensqr (Linear Algebra)
function comment givensqr (A,B)
  QR decomposition of A and B
  
  Returns {R=GA,Y=GB,c}, such that R is upper triangle, and c shows
  the linear independent columns of A.
  
  See: 
qrdecomp (Linear Algebra)
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.
  
  >A=normal(5,5);
  >{R,Q}=qrdecomp(A);
  >longest totalmax(abs(Q.A-R))
  4.440892098500626e-016
  
  see: givensqr, givensrot
function comment orthogonal (A)
  Compute an orthogonal basis of the columns of A
  
  This functions uses the Gram-Schmid method to compute an orthogonal
  basis of the columns of A. c is a row vector with 1 for the indices
  of the basis columns of A.
  
  >shortformat; A=random(3,5)
  0.883227  0.270906  0.704419  0.217693  0.445363
  0.308411  0.914541  0.193585  0.463387  0.095153
  0.595017  0.431184   0.72868  0.465164  0.323032
  >{B,c}=orthogonal(A); c
  [1,  1,  1,  0,  0]
  >O=B[,nonzeros(c)]
  0.796621 -0.370762 -0.477421
  0.278169   0.92606  -0.25502
  0.536672 0.0703505  0.840853
  >O.O'
  1         0         0
  0         1         0
  0         0         1
  
  See: 
svdimage (Linear Algebra)

Fits and Regression Analysis

There are functions for linear regression using the normal equation, orthogonal transformations, and singular value decomposition.

function fitnormal (A:complex, b:complex)
  Computes x such that ||A.x-b|| 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 (Numerical Algorithms)
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 minimal.
  
  The fit uses the normal equation solved with an orthogonal
  decomposition.
  
  w : Optional weights (squared) for the fit.
  

Norms

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)
  Maximal row sum of A
function comment scalp (v,w)
  Scalar product of v an w
function crossproduct (v:real, w:real)
  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. It will also work for a matrix of 1x3 columns
  v, and a matrix of 1x3 columns w.
  
  See: 
vectorproduct (Linear Algebra),
scalp (Euler Core),
scalp (Linear Algebra),
scalp (Geometry in Euler)
function vectorproduct (A:real)
  Multiplication of the columns of an nx(n-1) matrix A.
  
  Generalized cross product of n-1 vectors of length n. The result is
  perpendicular to the vectors, and its length is the volume of the
  polyhedron spanned by the vectors squared.
  
  A : real matrix of size n times (n-1)
  
  See: 
crossproduct (Euler Core),
crossproduct (Linear Algebra),
scalp (Euler Core),
scalp (Linear Algebra),
scalp (Geometry in Euler),
gramdet (Linear Algebra)
function gramdet (A)
  Gram determinant of the columns of A.
  
  The Gram determinant is the volume of the polyhedron spanned by the
  columns of A.
  
  A : real matrix of size kxn (k<n)
  
  See: 
vectorproduct (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 Gauss-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 (Linear Algebra)
function echelon (A)
  Forms A into echelon form using the Gauss algorithm
  
  See: 
LU (Linear Algebra),
lu (Linear Algebra)

Eigenvalues and Singular Values

For eigenvalues and eigenvectors, Euler has algorithms based on the AlgLib library, and algorithms based on the characteristic polynomial.

function comment charpoly (A)
  Characteristic polynomial of A
  
  Euler uses an orthogonal transformation to a simpler matrix, and a
  recursion to compute the characteristic polynomial. Note that
  polynomials in Euler start with the constant coefficient.
function eigenvalues (A:complex,
    usecharpoly=0, check=1)
  Eigenvalues of A.
  
  Returns a complex vector of eigenvalues of A.
  
  For real matrices, the functions uses an algorithm from AlgLib and
  LAPACK to find the eigenvalues, unless usecharpoly is true.
  
  If usecharpoly is true or the matrix is complex, 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 (Linear Algebra),
charpoly (Maxima Documentation),
svd (Euler Core),
jacobi (Linear Algebra),
jacobi (Maxima Documentation),
svdeigenvalues (Linear Algebra),
mxmeigenvalues (Maxima Functions for Euler),
aleigen (Linear Algebra)
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,
    usekernel=0, usecharpoly=0, check=1)
  Eigenvectors and a basis of eigenvectors of A.
  
  Returns {l,x}, where l is a vector of eigenvalues, x is a basis
  of eigenvectors.
  
  For real matrices A an algorithm from AlgLib and LAPACK is used,
  unless usekernel is true.
  
  usekernel : Compute the eigenvalues, then the kernels.
  usecharpoly : Use the characteristic polynomial for the
  eigenvalues.
  check : Check the return code of AlgLib.
  
  See: 
eigenvalues (Linear Algebra),
eigenvalues (Maxima Documentation),
eigenspace (Linear Algebra),
mxmeigenvectors (Maxima Functions for Euler)
function comment aleigen (A,flag)
  AlgLib routine to compute the eigenvalues of A
  
  This function returns the eigenvalues {v,res}, of the eigenvalues
  and the eigenvectors {v,res,W}, if flag is true. A must be real
  vector. The res flag is true, if the algorithm succeeds (from
  AlgLib and LAPACK).
  
  It is easier to use the functions eigenvalues() and eigen().
  
  See: 
eigenvalues (Linear Algebra),
eigenvalues (Maxima Documentation)
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 (Maxima Documentation)
function overwrite jacobi (A:real, eps=none)
  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)
function choleskyeigen (A, eps=none)
  Iterates the Cholesky-iteration, until the diagonal converges.
  
  A must be positive definite symmetric and at least 2x2.
  
  See: 
cholesky (Maxima Documentation)

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.
  
  See: 
svd (Euler Core)
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.
  This can be used to compute an orthogonal basis of vectors.
  
  >shortformat; A=random(3,2)
  0.493453  0.601344
  0.659461  0.967468
  0.193151  0.935921
  >svdimage(A)
  -0.441518  -0.45828
  -0.357159  -0.69891
  0.823103 -0.549094
  
  See: 
orthogonal (Linear Algebra),
svd (Euler Core),
kernel (Linear Algebra),
image (Maxima Documentation),
svdkernel (Linear Algebra)
function svdcondition (A:real)
  Condition number based on a singular value decomposition of A
  
  Returns the quotient 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
  
  See: 
svd (Euler Core)
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 eigenvalue 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

Gauss-Jordan Scheme

Some functions for the demonstration of the Gauss-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 Gauss-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
  
  See: 
gjstep (Linear Algebra)
function gjstep (A:numerical, i:index, j:index,
    n:positive integer vector,
    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_Jordan algorithm.
  
  The function modifies A and n.
  
  n : A row vector of indices, the function assumes the Gauss-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 making 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].
  
  You can add a vector v of strings, which contains the names of the
  variables. The problem will then be printed using gjprint.
  
  See: 
pivotize (Linear Algebra),
gjprint (Linear Algebra)
function gjsolution (M : real, n : positive integer vector)
  Read the values of the solution from the Gauss-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 Gauss-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}

Additional Matrix Functions

function comment magic (n)
  Magic square of size nxn.

Documentation Homepage