8.   Linear Algebra

2015-05-04

$Version

8.1  Constructing Matrices

M = list = {{a11,a12,a13,...a1n}, {a21,a22,a23, ..., a2n}, 
            ....,                 {am1,am2,am3, ..., amn} }
            
M = Table[f[i,j], {i,m}, {j,n}]		build a  (n x m) matrix, whose element 
							M[i,j] is given by the value f[i,j]
					
A = Array[a, {m,n} }  	build a  (n x m) matrix, whose element a[i,j]
						  is the symbol  a[i,j]

DiagonalMatrix[list]		generate a diagonal matrix with the elements of   
					list  on the main diagonal
                             
IdentityMatrix[n]   		generate the  (n x n)  Identity Matrix

TableForm[list]	      		Prints  list  with the elements of  list  arranged  in a regular array

MatrixForm[list]	   		  Prints  list  as a matrix

A = {{1,2,3,4}, {2,3,0,-5}, {2,-1,1,1},{-2,2,0,-5} }

MatrixForm[A]

M = Table[m[i,j], {i,2}, {j,3}]

MatrixForm[M]

f[i_,j_] = i/j;

MF = Table[f[i,j], {i,3}, {j,3} ]

TableForm[MF]

Clear[a,aa]
aa = Array[a,{3,4}]

MatrixForm[%]

a[i_,j_]:=ToExpression["a"<>{ToString[i] ,ToString[j]}]

MatrixForm[aa]

Clear[a,aa]
aa = Array[a,{3,4}];

a[i_,j_] := Subscript[a , i,j]

MatrixForm[aa]

Di = DiagonalMatrix[{2,1,0,-1,-2}];
MatrixForm[Di]

MatrixForm[IdentityMatrix[3]]

MatrixForm[0 IdentityMatrix[3]]

MatrixForm[Table[0, {2}, {4}]]

Array[0, {2,4}]

MatrixForm[Table[If[i > j,0,1], {i,3},{j,3}]]

MatrixForm[Table[If[i == j,0,1], {i,3},{j,3}]]

MatrixForm[Table[If[i <= j,0,1], {i,3},{j,3}]]

8.2  Getting Pieces  of  Matrices

m[[i,j]]			the  i,j -th  element

m[[i]]  			the  i -th  row 
  
Transpose[m][[i]]  	the  i -th  column
  
m[[{i1,...,ir},{j1, ..., js}]]	    the  r x s  submatrix of  m  with elements having 
						          row  indices  ik   and column indices  jk

m[[ Range[i0,i1], Range[j0,j1] ]]     a submatrix with elements  having row and column
						           indices respectively,  in the ranges given
						         		
Diagonal[m]		gives a list containing the elements of the leading diagonal of the matrix  m						         	
Diagonal[m, k]		gives a list containing the elements of the  k-th  diagonal of the matrix  m

Diagonal works also on rectangular matrices.

A = {{1,2,3,4}, {2,3,0,-5}, {2,-1,1,1},{-2,2,0,-5} }; MatrixForm[A]

Diagonal[A,-3]

A[[3,2]]

A[[2]]

B = Transpose[A]; MatrixForm[B]

B[[2]]

MatrixForm[ A[[ {2,4}, {1,3} ]]  ]

MatrixForm[  A[[ Range[2,4], Range[1,3] ]]  ]

K = Drop[A,1];
MatrixForm[K]

Diagonal[K]

Diagonal[K,1]

8.3   Scalars, Vectors, Matrices

Mathematica  represents matrices and vectors using lists. Anything  that is not a list is considered  a scalar.  
A vector  consists of a list of scalars. A matrix consists of a list of vectors, representing  each of its rows; each of these vector s must have the same length.
Operators for testing the structure of vectors and matrices:

VectorQ[expr]		= True    	expr  has the form of a vector
               		= False    	otherwise
                
MatrixQ[expr]  	 = True      	expr  has the form of a matrix
              		= False 	otherwise
                              
Dimensions[expr]  			a list of the dimensions of a vector  or of a matrix

v = {a, b, c} ;  VectorQ[ v ]

VectorQ[ p + q  ]

Di = DiagonalMatrix[{2,1,0,-1,-2}];

Dimensions[Di]

Dimensions[v]

MatrixQ[ {Di, v} ]

MatrixQ[ Di ]

8.4  Operations on Scalars, Vectors, Matrices

Log[ {a, b, c} ]

MatrixForm[ Exp[Di] ]

D[ {1, x, x^2, x^3}, x]

{a, b} + {c, d}

{a, b, c }  + {d, s}

Thread::"tdlen" : Objects of unequal length in {a, b, c }  + {d, s}cannot be combined. >>

In general, addition and subtraction can be applied to like lists only.  However, a scalar is interpreted as unit list having the same dimensions as the list representing the other operand.

1 + {a, b}

{a, b} + p

% /. p -> {c, d}

Exp[Di] + c; MatrixForm[%]

k {a, b}

{a, b}/k

A//MatrixForm

Diagonal[A]

Diagonal[A,2]

Diagonal[A,-1]

D = Drop[A,1];
MatrixForm[%]

Diagonal[%% ]

Diagonal[%%%,1 ]

8.5  Multiplying  Vectors and Matrices

Different kinds of vector  and matrix multiplication:

s  =  scalar;    v  =  vector;    m  =  matrix.

s v, s m   			multiply each element by a scalar         
                           
v.v, v.m, m.v, m.m     	vector and matrix multiplication

v = {a, b, c }

s v

v . {ap, bp, cp}

v . v

{a,b,c} {d, e, f}

{{a, b}, {c, d}} . {x, y}

m1 = {{a, b}, {c, d} }; m2 = {{1,2}, {3,4}};

MatrixForm[m1]

MatrixForm[m2]

m = m1 . m2

MatrixForm[m]

v = {x, y}

MatrixForm[ m1 . v ]

MatrixForm[ v . m1]

ExpandAll[ MatrixForm[ v . m1 . v] ]

8.5.1 The  outer product

Outer[Times, t1, t2, ...]    	Outer Product

Outer[Times, {a, b, c}, {x, y, z} ]

MatrixForm[%]

The command  Outer[Times,v1,v2]  does not give the vector product of two vectors  but the tensor  product.

r = {x, y, z};

IdentityMatrix[3] r.r - Outer[Times, r, r]//MatrixForm

Moment of inertia = second moment =  Trägheitstensor

The operator  in Outer[] need  not be Times. It  may be any function f. The lengths  of the  lists  need  not be equal.

f[i_,j_]=i/j;
Outer[f,{1,2,3},{2,4,8,10}]

Outer[List,{1,2,3},{2,4,8,10}]

Outer[Times, {a, b}, {c, d}, {e,f} ]

MatrixForm[%]

Clear[x,y,z];    r = {x, y, z};

m1 = {{a, b}, {c, d} }; m2 = {{1,2}, {3,4}};

m12 = MatrixForm[ Outer[Times, m1, m2] ]

m21 = MatrixForm[ Outer[Times, m2, m1] ]

8.5.2  The Kronecker product

The Direct Product (or Kronecker Product) of two matrices is defined as (where  there are  no  submatrices !)

A  ⊗   B      =     (({{(Subscript[a, 1,1])B, (Subscript[a, 1,2])B     ...,  (Subscript[a, 1,n ])B}, { (Subscript[a, 2,1])B,  (Subscript[a, 2,2]) B    ...,                (Subscript[a, 2,n])B
		....}, { (Subscript[a, m,1])B,  (Subscript[a, m,2])  B  ...,  (Subscript[a, m,n])B}}))

So the Outer Product of Mathematica does not correspond to the Kronecker  product. However, there is a special command KroneckerProduct[] :

Clear[m1,m2]
m1= {{a,b},{c,d}};
m2 = {{1,2,3},{4,5,6}};
MatrixForm[m1]
MatrixForm[m2]

KroneckerProduct[m1,m2]//MatrixForm

8.6 Basic Matrix Operations

Transpose[m]         transpose of   m
  
Inverse[m]           inverse of  m
 
Det[m]               	 determinant of  m

Minors[m, k]          a matrix of the  k x k  minors of  m

MatrixRank[m]	rank  of the matrix m 

Tr[m]     		  trace  of  m

MatrixForm[A]

MatrixForm[ Transpose[A] ]

Det[A]

B = Inverse[A]; MatrixForm[B]

MatrixForm[ A . B ]

G = A; G[[1]] = G[[2]]; MatrixForm[G]

So  G  is a singular  matrix:

Inverse[G]

Inverse::sing : Matrix {{2, 3, 0, -5}, {2, 3, 0, -5}, {2, -1, 1, 1}, {-2, 2, 0, -5}}  is singular. >>

Det[G]

MatrixRank[G]

ma = Array[a,{3,4}]; MatrixForm[ma]

mi3 = Minors[ma,3]

Dimensions[mi3]

mi2 = Minors[ma, 2]

Dimensions[mi2]

mi1 = Minors[ma,1]//MatrixForm

Minors[G,3]

So the matrix  G  has  rank  3.

Tr[A]

Tr[G]

8.6.1  Characteristic Polynomial of a Matrix   A :

A = {{1,2,3,4}, {2,3,0,-5}, {2,-1,1,1},{-2,2,0,-5} };

f=Det[A-x IdentityMatrix[4]]

cp = CharacteristicPolynomial[A, λ]

Solve[f == 0., x]

8.6.2  Functions of a Matrix  A

MatrixPower[m, n]    	n -th  power of matrix  m
  
MatrixExp[m]          	matrix exponential

Clear[a,b,c,d]
m = {{a,b}, {c,d}}; MatrixForm[m]

MatrixForm[ m^2 ]

MatrixForm[ m . m ]

MatrixForm[ MatrixPower[m, 2] ]

dd = DiagonalMatrix[{1,-1}]; MatrixForm[dd]

MatrixExp[dd]//MatrixForm

dd = I x DiagonalMatrix[{1,-1}];MatrixForm[dd]

MatrixExp[dd]

ComplexExpand[%]// MatrixForm

MatrixExp[m]

FullSimplify[%]

8.7 Solving Linear Systems

Systems  of linear equations can be solved by various Mathematica  commands.

m  =  matrix,  x  =  vector (unknowns),  b  =  vector (right hand side of an equation) .

						m.x  == b

Solve[eq, unk]	     Solves the system of linear equations   eq  for the unknowns contained  in
			         list   unk. Works for regular and singular systems.

LinearSolve[m, b]      gives a vector x  which solves the matrix equation  m.x == b 
                     
NullSpace[m]    a list of basis vectors whose linear combinations x  fulfil the matrix  eqation  m.x == 0                      

RowReduce[m]     a simplified form of m  obtained by making linear combinations of rows

Let  n  denote the number of unknowns = number of columns of the matrix m. 
r  is the rank of the matrix (command MatrixRank[], § 8.6 ). 
If  m is a square matrix and is regular, n = r, then the system of equations has a unique solution, in particular  x = 0  if b = 0.
If r < n the homogeneous system does not have a unique solution; in fact, there is a solution manifold depending on  n - r  arbitrary parameters. The command NullSpace[] gives a basis for that manifold. 
If r > n, there is no exact solution; an approximate solution may be found by optimization methods.
The inhomogeneous system with r < n has solutions only if the vector b on the right hand side does not introduce a contradiction.
This implies: the rank of the matrix formed by adjoining b to m is equal to the rank of m; this condition is equivalent to the following one: all vectors of the nullspace of the matrix (OverHat[m]) (= matrix transposed to m) are orthogonal to b.

m = {{1, 5}, {2, 1}}; MatrixForm[m]

The commands  listed above provide several  methods for solving systems  of linear equations.

8.7.1	Method 1    uses the command   Solve[].

m . {x, y} == {a, b}

Solve[%, {x,y}]

Solve[]    works also for solvable singular linear systems:

{x + y, 2x + 2y} == {0,0}; 
Solve[%, {x,y}]

Solve::"svars" :  Equations may not give solutions for all "solve" variables. >>

The following singular system is contradictory;  thus no solution exists:

{x + y, 2x + 2y} == {a,b}; 
Solve[%, {x,y}]

8.7.2	Method  2  for solving a regular system uses the inverse  matrix :

Inverse[m].{a,b}

8.7.3	In  method  3   one works with the command    LinearSolve[].

LinearSolve[m, {a,b} ]

LinearSolve[] does not give the general solution of a solvable  singular system, but only a special solution.

LinearSolve[{{1,1},{2,2}}, {0,0} ]

For  large systems  methods  3 and 4 work faster  than  the first two methods.

8.7.4	In method 4  one uses the command   RowReduce[].

For  that  purpose  the vector  on the right hand side must be joined  to the matrix of coeffcients  to give a 
rectangular  matrix:

h1 = Append[m[[1]], a];
h2 = Append[m[[2]], b];
ma = {h1,h2}; MatrixForm[ma]

RowReduce[ma] //MatrixForm

The last column of the above matrix contains the solutions of the linear  system.

mm = {{1,1,0}, {2,2,0}}; MatrixForm[mm]

RowReduce[mm]//MatrixForm

mm = {{1,1,a}, {2,2,b}}; MatrixForm[mm]

RowReduce[mm]//MatrixForm

This shows again that the singular system with the rhs {a,b} is contradictory.

Method  4   is shown again for the following system:

			- 3 x  +  2 y  + 11 z   =   1,
			      x  +  3 y  +   7 z  =   -5,
			  -2 x  -  3 y   +   5 z  =   2.

g = {{-3,2,11,1},{1,3,7,-5},{-2,-3,5,2}}; MatrixForm[g]

rr = RowReduce[g]; MatrixForm[rr]

v = Transpose[rr][[4]]

This is the solution; as is verified below:

cm = g[[ Range[3], Range[3] ]]; MatrixForm[cm]

cm . v == Transpose[g][[4]]

me = {{1,1,1,-1}, {1,2,3,-4}, {1,3,6,-10},{1,4,10,-a}}

me[[Range[4],Range[3]]] . {x,y,z} == me[[4]]//Thread

Solve[%, {x,y,z}]

me[[ Range[3],Range[3]] ] . {x,y,z} == me[[ 4,Range[3] ]]//Thread

Solve[%, {x,y,z}] //Flatten

So the system  of  4  equations  above  is overdetermined.  But the system  comprising  the first  3  equations 
is solvable.  This is also verified  using   RowReduce[].

met = RowReduce[me]; MatrixForm[met]

mer = RowReduce[ me[[Range[3],Range[3] ]] ]; MatrixForm[mer]

Now we look for the value of   a   for which the system of   4   equations is  solvable:

sa = Solve[Det[me] == 0, a]//Flatten

RowReduce[me /. sa]//MatrixForm

Solve[(me /. sa)  . {x,y,z,1} == 0, {x,y,z}]

8.7.5	Eliminate

??Eliminate

eqs={2x+3y+z==1,x-y-z==4,3x+7z==5};

Solve[eqs,{x,y,z}]//Flatten

el=Eliminate[eqs,z]

Apply[List,el]

Solve[%,{x,y}]//Flatten

8.7.6	The Null Space

A regular matrix has no nullspace.

A = {{1,2,3,4}, {2,3,0,-5}, {2,-1,1,1},{-2,2,0,-5} };
MatrixForm[A]

Det[A]

NullSpace[A]

Now a singular matrix is constructed:

G = A;  G[[4]] = G[[1]]; MatrixForm[G]

rg = MatrixRank[G]

Length[G] - rg

v = NullSpace[G]//Flatten

G . v

v . G

v.Transpose[G]

x = {x1, x2, x3, x4}

so = Solve[G.x == 0, x]//Flatten

Solve::svars : Equations may not give solutions for all"solve"variables. >>

xs = x /.so /. x4 -> 25  α

v  α

G.xs

A new singular matrix is constructed:

H = G;H[[3]] = 2 G[[1]]; MatrixForm[H]

rh = MatrixRank[H]

Length[H] - rh

so = Solve[H.x == 0, x]

Solve::svars : Equations may not give solutions for all"solve"variables. >>

ns = NullSpace[H]

ys = α ns[[1]] + β ns[[2]]

H.ys//Together

Now a soluble inhomogeneous singular system H.x = b  is constructed.

b = {b1,b2,b3,b4}

This vector b  must be orthogonal to all solutions of the homogeneous system transposed to H:

nt = NullSpace[Transpose[H]]

sys = {nt[[1]]. b, nt[[2]]. b}

so = Solve[sys == 0,b] //Flatten

vb = b/. so

Only such vectors on the right hand side of   H.x = b  lead to a non-trivial solution:

Hb = Append[Transpose[H],vb /. so]//Transpose; MatrixForm[Hb]

MatrixRank[H]
MatrixRank[Hb]

8.8  Eigenvalues and Eigenvectors. Jordan Decomposition

An eigenvector  X  of a square   n x n  matrix  A  is a vector  whom the application of the matrix transforms  into a collinear vector  α X such  that

						AX = α X .					(1)

Any vector fulfilling this condition is called an eigenvector   (Subscript[X, i]) ;   the corresponding  value  α  is called the eigenvalue  (Subscript[α, i]) .
The above equation is rewritten by use of the identiy matrix E:

					(A  -  α E)X  =  0.					(2)
						
This linear system for the components of the vector   X  has a nontrivial solution only if the matrix of the coefficients is singular. The condition that its determinant must be zero, gives a polynomial of order  n  in  α, the characteristic  polynomial  P(α) , the characterstic equation P(α) = 0 respectively:

P(α) := det(A  -  α E)  =  (-1)^n[ (α^n) -  ((Subscript[a, 11])+  (Subscript[a, 22]) + ... + (Subscript[a, 1nn])) (α^(n-1)) + - ... +  (-1)^ndet(A)]  =  0.	(3)

Its zeros are the eigenvalues   (Subscript[α, i]).  After an eigenvalue  (Subscript[α, i])  has been inserted for α  in eq.(2) one can solve it to find the corresponding eigenvector   (Subscript[X, i]).

A matrix  M  may be transformed  by an equivalence  transformation.  Such a transformation  is performed with the help of a regular matrix  T such that 

					B  =  (T^(-1)) A T.
						
If the matrix  T  consists  of the eigenvectors  (Subscript[X, i]) of the matrix  A, then the matrix B becomes a diagonal matrix, i.e.
					(T^(-1)) A T  =  diag((Subscript[α, 1]), (Subscript[α, 2]), ..., (Subscript[α, n])) .

Multiple eigenvalues are called degenerate. Then it may happen that  the number of independent eigenvectors is smaller than the dimension of the matrix A. In this case it is impossible to build the matrix  T  and to transform the Matrix A into a purely diagonal form. But, it can be transformed to Jordan normal form (cf. next section).

The following types of matrices always can be diagonalized:
1. Matrices whose eigenvalues are pairwise different, i.e. all the eigenvalues are non-degenerate.
2. Normal matrices; they fulfil the following condition: A(SuperPlus[A]) =  (SuperPlus[A])A, where (SuperPlus[A])= (SuperStar[(OverTilde[A])]) is the complex conjugatetranspose to A.
2.1  Real (A = (SuperStar[A])) normal matrices are the following ones:
	2.1.1 symmetric (A = (OverTilde[A])), 2.1.2 antisymmetric (A = - (OverTilde[A])), 2.1.3 orthogonal (A(OverTilde[A]) = E), 
	2.1.4 anti-orthogonal (A(OverTilde[A]) = E);
2.2  Complex normal matrices are the following ones:
	2.2.1 Hermitian (A =  (SuperPlus[A])), 2.2.2Anti-Hermitian (A = - (SuperPlus[A])), 2.2.3 unitary  (A(SuperPlus[A])= E),  
	2.2.4 antiunitary  (A(SuperPlus[A])= -E), 
The eigenvectors of real matrices of type 2.1.1 and 2.1.2 can be chosen to be real and  orthonormal; then the matrix  T  is orthogonal. In all other cases the eigenvectors  can be chosen to be complex  orthonormal; then the matrix  T  is unitary.

8.8.1 Eigenvalue Decomposition of   a  Matrix

Eigenvalues[m]		a list of the eigenvalues  of  m

Eigenvectors[m]	   	a list of the eigenvectors  of  m

Eigensystem[m]       	 a list of the form {eigenvalues,  -vectors}
                                       
Eigenvalues[ N[m]]    	numerical eigenvalues  

Eigenvalues[ N[m,k]]  	numerical eigenvalues,  starting  with k -digit precision

Clear[x]

A = {{1,1,2}, {1,2,1}, {2,1,1}}; MatrixForm[A]

Eigenvalues[A]

The characteristic  equation  ceq  of the matrix is is found below:

AM = A - x IdentityMatrix[3]; MatrixForm[AM]

ceq = Det[AM]

ds = Solve[ceq == 0,x]

Eigenvectors[A]

Eigensystem[A]

eigenvalues       	vector1      	vector2    	 vector3

AM1 = AM /. ds[[1]]

NullSpace[AM1]

AM /. ds

NullSpace[ %[[2]] ]

Clear[a,b,c,d,m]
m = {{a,b},{c,d}}; MatrixForm[m]

Eigenvalues[m]

Eigensystem[m]

8.8.2 Tranforming a Matrix to Jordan Form

Not every  matrix can be diagonlized. see the discussion preceeding the previous subsection.   Matrices  with multiple eigenvalues may have less eigenvectors  than eigenvalues. These may be transformed  to a normal form with 1's in the diagonal above the main diagonal containing the degenerate eigenvalues.  The JordanDecomposition[] applied to a diagonalizable  matrix yields the same diagonal matrix as the eigenvalue decomposition.

JordanDecomposition[m]        yields the Jordan decomposition of a square matrix  m.  The result is a list{s, j}    where  s  is a similarity  matrix  and  j   is the Jordan canonical  form of m.

mm = {{1,0,0},{2,1,0},{0,0,-1}}; MatrixForm[mm]

Eigensystem[mm]

Note the last vector  in the list above is a Null vector.

jd = JordanDecomposition[mm]

MatrixForm[jd[[2]]]

J = Inverse[jd[[1]]].mm.jd[[1]] ; MatrixForm[J]

8.9  Singular Value Decomposition.   PseudoInverse

In this paragraph generalizations  are treated  which are analogous  to eigenvalue decomposition and inverse matrices, but for rectangular  and even singular matrices.
The singular value decomposition is applicable to rectangular matrices; it is a procedure which somewhat (but not exactly) resembles the diagonalization (=  eigendecomposition )  of a square matrix. In fact:
Any real m x n  (with m ≥ n)  matrix   A  has a representation as  A = (Overscript[u, ∼]) w v  where  (Overscript[u, ∼])u = E, (Overscript[v, ∼])v = E and w is a diagonal matrix with non-negative entries. The matrix  u  is  m x m, v is  n  x n  and  w  is  m x n. The singular values (Subscript[w, i]) are the square roots of the eigenvalues of  (Overscript[A, ∼])A, and the columns of  v   are the eigenvectors of  (Overscript[A, ∼])A. The colums of  u  are the eigenvectors  of  A(Overscript[A, ∼]).
Computation of the singular value decomposition is more stable than forming   (Overscript[A, ∼])A   and computing the eigendecomposition.

8.9.1 Singular Value Decomposition

SingularValueDecomposition[m]
SingularValues[m] gives the singular  value decomposition  for a numerical matrix m. 
			 The result is a list {u, w, v}, where w is the list of singular values, and  
	m  can be written as     m = u. w.  Conjugate[Transpose[v]]  ;     or as   
	                                     m = Conjugate[Transpose[u]].w.v .
 	
 SingularValueList[m]    gives a list containing the singular values of the matrix   m.

This method is applicable to rectangular matrices, too. Even if the matrix  m  is square and diagoni-zable, the singular values  contained  in  w  differ from the eigenvalues.    
If   u ( v )   is square,  it is orthogonal  (or unitary if the original matrix  m  is complex) ; otherwise it is orthogonal (or unitary) in one direction only.

Example 1

A = {{1,1,2}, {1,2,1}, {2,1,1}}; MatrixForm[A]

SingularValueList[A]

The singular values contained in the diagonal matrix  w  shown in the above list  differ from the eigenvalues of  A  given below:

Eigenvalues[A]

{u,w,v}=SingularValueDecomposition[A]

MatrixForm[u]

MatrixForm[w]

MatrixForm[v]

u.w.Transpose[v]//MatrixForm

This is the original matix.

Transpose[u].w.v//MatrixForm

Example 2

A = {{1,1,2}, {1,I,1}, {2,1,1}}; MatrixForm[A]

{u,w,v}=SingularValueDecomposition[A]

w

N[w]

Eigenvalues[A]

N[%]

A == u.w.Conjugate[Transpose[v]]//Simplify

Example 3: The singular value decomposition of a rectangular matrix

ma={{1,-2},{2,-1},{1,1}};MatrixForm[ma]

{u,w,v} = SingularValueDecomposition[ma]//Chop

MatrixForm[w]

u.w.Transpose[v] == ma

Example 4

ma = {{1,0,0},{2,1,0},{0,0,-1}};MatrixForm[ma]

{u,w,v}=SingularValueDecomposition[ma];

MatrixForm[u]

MatrixForm[w]

MatrixForm[v]

MatrixForm[u.w.Transpose[v]]//FullSimplify

In this particular  case  the elements  of all three  matrices  may be expressed   by symbolic expressions.

wa= {(Sqrt[2]) + 1,1,(Sqrt[2])-1}; 
N[wa] - w //Chop

Eigenvalues[ma]

Compare these  eigenvalues  to the elemnts of  w  or  wa. The elements of the transformation matrix are:

v1 = (Sqrt[(1/2)-(1/(2 (Sqrt[2])))])               (* =Sin[π/8] *); 
v2 = (Sqrt[(1/2)+(1/(2 (Sqrt[2])))])               (* = Cos[π/8] *);

ua = {{v1,0,-v2},{v2,0,v1},{0,-1,0}};MatrixForm[ua]

ua//N
u//N

va = {{v2,0,-v1},{v1,0,v2},{0,1,0}}; MatrixForm[va]

va//N
v//N

Example 5

ma = {{1,0,0,3},{2,1,0,1},{0,0,-1,1}};MatrixForm[ma]

{u,w,v}=SingularValueDecomposition[ma];

MatrixForm[w]//N

MatrixForm[u]//N
MatrixForm[v]//N

MatrixForm[u.w.Transpose[v] - ma//N//Chop]

Transpose[u].u//N//Chop//MatrixForm

u.Transpose[u]//N//Chop//MatrixForm

v.Transpose[v]//N//Chop//MatrixForm

Transpose[v].v//N//Chop//MatrixForm

8.9.2 Pseudoinverse

The inverse (A^(-1)) of a matrix A exists  only if A is square  and has full rank. In this case, Ax = b  has the solution  x = (A^(-1))b. The pseudoinverse   (OverBar[A])  is a generalization  of the inverse, and  exists for any m × n matrix. We  assume   m > n.  If A has full rank (n) we define:

					 (OverBar[A]) =  ((OverTilde[A]) A)^(-1) (OverTilde[A])				    		(1)
						
with   (OverTilde[A])  the  transpose to  A, and the solution of    Ax = b is 

						x = (OverBar[A]) b. 					(2)

The best  way to compute  (OverBar[A]) is to use singular value decomposition. 
If W = ((Subscript[w, ij ])) is a diagonal matrix with the entries  wi, then we define the pseudoinverse  (OverBar[W]) as:

	                                                               (OverBar[W])   =  ((OverBar[(Subscript[w, ii])]))	=  (1/wi)     if   wi ≠  0 ,
						=   0       otherwise.
	
With, where U  (m×m) and V (n×n)  are orthogonal and  W(m,n) is diagonal with real, non-negative singular values wi, we get the pseudoinverse (OverBar[A])  to A  as:

					A =  UW(OverTilde[V])	→      (OverBar[A])  =   V(OverBar[W])(OverTilde[U])

For  nonsingular square  matrices  the pseudoinverse  agrees  with the common inverse.  But the pseudoinverse exists  also for singular  matrices.

If the rank r of A is smaller  than  n, r < n, then the inverse of A does not exist, and  one uses  only the first  r  singular values;  W   then becomes  an (r x r) matrix and U,V shrink accordingly.

Sources: Rudolf K. Bock, Pseudoinverse.  http://rkb.home.cern.ch/rkb/AN16pp/node220.html
	   http://www.das-gelbe-rechenbuch.de/download/Swz.pdf

PseudoInverse[m]	compute the pseudoinverse   of the rectangular  matrix  m; 
				works for numeric and symbolic matrices,  real or complex.

Example 1

ma={{1,-2},{2,-1},{1,1}};
mb={-1,1,5};
lx={x,y};
eq=ma.lx==mb//Thread

This is an overdetermined system; the coefficient matrix is rectangular, viz.  3 x 2. Its rank is 2.

MatrixRank[ma]

Eqs. (1)  and (2) are now illustrated:

pa = Transpose[ma].ma

ba = Inverse[pa]. Transpose[ma]

ba.mb

pia=PseudoInverse[ma]
pia.ma //MatrixForm
pso=pia.mb

Example 2

ma = {{1,2,3},{I, 1. ,-I}}//Transpose; MatrixForm[ma]

pa = PseudoInverse[ma]; MatrixForm[pa]

pa.ma//Chop//MatrixForm

Inverse[Conjugate[Transpose[ma]].ma] .Transpose[ma]

Example 3

Clear[a,b,c,d,m]
m = {{a,b},{c,d}}; MatrixForm[m]

Inverse[m]//MatrixForm

PseudoInverse[m]//FullSimplify//MatrixForm

MatrixForm[m = {{1,2},{1,2}}]

Inverse[m]

Inverse::sing : Matrix {{1, 2}, {1, 2}} is singular . >>

MatrixForm[p = PseudoInverse[m]]

This  shows that  p  satisfies the four conditions  that  uniquely characterize  the pseudoinverse p   of a matrix m :

m.p.m == m
p.m.p == p
p.m == Transpose[p.m]
m.p== Transpose[m.p]

mm= {{1,2,3},{3,2,1}};

MatrixForm[p = PseudoInverse[mm]]

Geometric  interpretation of the solution by the pseudoinverse

posize=PointSize[0.02];

3 independent  straight  lines

c1={1,1};
c2={3,2};
c3={2,3};
cp=Map[Point,{c1,c2,c3}]

g1=x-2y==-1
g2=2 x-y==1
g3=x+y==5

s1 = y /.Solve[g1,y]//Flatten
s2 =y /.  Solve[g2,y]//Flatten
s3 = y /. Solve[g3,y]//Flatten

gr1 = Plot[{s1,s2,s3}//Flatten,{x,0,3.2},PlotRange->{0,3.2},Ticks->{Range[3],Range[3]},Epilog->Prepend[cp,posize]]

ma={{1,-2},{2,-1},{1,1}};
mb={-1,1,5};
lx={x,y};
eq=ma.lx==mb//Thread

pia=PseudoInverse[ma]
pia.ma //MatrixForm
pso=pia.mb

Show[gr1,Graphics[{Hue[0],posize,Point[pso]}]]

The "pseudosolution" is the centre  of gravity of the triangle.

2 parallel straight lines crossing  a third one

c1={1,1};
c2={3,2};
c3={2,3};
cp=Map[Point,{c2,c3}];

g1=x-2y==-1;
g2=x-2 y==-3.95;
g3=x+y==5;

s1 = y /.Solve[g1,y]//Flatten
s2 =y /.  Solve[g2,y]//Flatten
s3 = y /. Solve[g3,y]//Flatten

gr2 = Plot[{s1,s2,s3}//Flatten,{x,0,3.2},PlotRange->{0,3.2},Ticks->{Range[3],Range[3]},Epilog->Prepend[cp,posize]]

ma={{1,-2},{1,-2},{1,1}}; MatrixForm[ma]
mb={-1,-3.95,5};
eq=ma.lx==mb//Thread

pia=PseudoInverse[ma]; MatrixForm[pia]
pia.ma//MatrixForm
pso=pia.mb

Show[gr2,Graphics[{Hue[0],posize,Point[pso]}]]

The "pseudosolution" is the centre  of gravity of the two crossing points..

3 parallel  straight  lines

g1=x-2y==-1;
g2=x-2 y==-3.95;
g3=x-2 y==-2.1;

s1 = y /.Solve[g1,y]//Flatten
s2 =y /.  Solve[g2,y]//Flatten
s3 = y /. Solve[g3,y]//Flatten

gr3 = Plot[{s1,s2,s3}//Flatten,{x,-1,3.2},PlotRange->{0,3.2},Ticks->{Range[-1,3],Range[3]},Epilog->Prepend[cp,posize]]

ma={{1,-2},{1,-2},{1,-2}};
mb={-1,-3.95,-2} ;
eq=ma.lx==mb//Thread

pia=PseudoInverse[ma]
pia.ma //MatrixForm
pso=pia.mb

Show[gr3,Graphics[{Hue[0],posize,Point[pso]}]]

8.10      Exercises

Only the ewxercises with a cross (+) are obligatory. If an obligatory exercise comprises several tasks distinghuished by a letter, then only one of these task is requested.

8.1+ 	Solve the following systems of linear equations;   3) and 4) are homogeneous systems, the matrix 
	gives the coefficients.  	Use at least 3 methods when the system is regular,  otherwise at  least  2 		methods.

a)	x + y + z = 1,	x + 2y + 3z = 4,    x + 3y + 6z = 10,	   x + 4y +   6z = 19.
	b)  	x + y + z = 1,    x + 2y + 3z = 4,    x + 3y + 6z = 10,       x + 4y + 10z = 19.

	c)        (({{0, 1, -3, -1}, {1, 0, 1, 1}, {3, 1, 0, 2}, {1, 1, -2, 0}}))		d)	  (({{1, 3, 2, 4}, {5, 2, 0, 1}, {3, -4, -4, -7}, {-7, 5, 6, 10}})) .

8.2+   	Replace  number  19   in  the system  8.1.b)  above by  an unknown   a. Without using the results 		obtained in  8.1.b)   find  the  value of   a   such  that the system of  4   equations has  
	a unique solution.

8.3+   	Compute the square, the transpose, the trace, the determinant, the inverse,  the characteristic 			polynomial, the latent roots (= eigenvalues) and eigenvectors of the following  matrices and  
	transform  them into  equivalent  diagonal  matrices. 

	a)       (({{10, -14, -10}, {-14, 7, -4}, {-10, -4, 19}})) ;   b)   (({{5, 2, 2}, {2, 2, 1}, {2, 1, 1}}));       c)   (({{-2, 0, -4}, {0, 2, 4}, {-4, 4, 0}}));    d)   (({{1, 0, (Sqrt[3])}, {0, 1, 1}, {(Sqrt[3]), 1, 1}})) ;

	e)        (({{1, -2, 1}, {-1, 2, -1}, {2, 1, 1}}));           f) *)     (({{8, 2, 3, 4}, {2, 5, 4, 5}, {3, 4, 5, 6}, {5, 6, 7, 9}}))

*) It is not reasonable to try to calculate analytic expressions for the eigenvectors of this matrix.

8.4+   	Some powers of the matrices given below coincide with the identity matrix; find the 
	correspondig lowest positive power.   

	a)     (({{-(1/2), ((Sqrt[3])/2)}, {((Sqrt[3])/2), (1/2)}}))              b)     (({{-(1/2), -((Sqrt[3])/2)}, {((Sqrt[3])/2), -(1/2)}}))            c)      (({{(1/(Sqrt[2])), -(1/(Sqrt[2]))}, {(1/(Sqrt[2])), (1/(Sqrt[2]))}}))


	d)        (({{0, 1, 0}, {0, 0, 1}, {-1, -1, -1}}))              e)   (({{0, 1, 0}, {0, 0, 1}, {1, 0, 0}}))

8.5   	a) Write a program, which tests all powers (up to a given limit) of a given matrix   a  to find 		the lowest positive power of   a  which equals the identity matrix;   
             b)  Write a program, which tests all powers (up to a given limit) of a given matrix   a  to find 	the lowest positive power of   a  which equals the null matrix  (then  a   is  nilpotent).

8.6  	Write a  program, which proves by calculation that a given matrix   a  satisfies its own 		characteristic equation.  (For this task commands described in Chap. 9  may be helpful ).

8.7+  	A zero divisor of a given matrix  A  is a (square or rectangular)  matrix L  or  R, which differs from 
	the null matrix, but for which   LA  =  0  (then L  is a left zero divisor)  or  AR = 0  (then  R  is a right 	zero divisor).  R, L  respectively,  consist of the null space of  A,  of the transpose of A  respectively.  	Find the left and right  zero divisors of the following matrices:

a)             (({{2, 0, -1}, {3, 4, 2}, {0, -8, -7}}))       ,

b)             (({{1, 1, 1, 1}, {1, 2, 3, 4}, {3, 5, 7, 8}, {6, 5, 4, 2}}))  .

8.8.  	Determine the shape (determined by the eigenvalues) and the orientation (given by the eigenvectors) 
	of the following conics:

	1)	   x 2 +   5 y 2   +   3 z 2  +  8 y z  -    8 x z 	= 1;	

	2)*	5 x 2 + 13 y2  + 14 z 2   -  4 x z  +  12 x y + 12 y z     	= 1.
	*) It is not reasonable to calculate the eigenvalues or eigenvectors as analytic expressions .

8.9.  	Compute the Kronecker product  of each of the 3 Pauli spin matrices and of the 			corresponding 	 2 x 2	identity matrix  e.

(Subscript[σ, x])  =    (({{0, 1}, {1, 0}})),         	σy  =   (({{0, -i}, {i, 0}})),         	σz =     (({{1, 0}, {0, -1}})).

8.10 	A matrix can be diagonalized if it is normal, i.e. if it  commutes with its adjoint ( A A+ = A+ A), 
	or has pairwise different eigenvalues. Write a program which checks this property for a given 	matrix A.  Apply this test to the following matrices:

a)          (({{1, 1}, {0, 2}}))

b)          (({{1, -i}, {i, 2}}))

8.12+  	Is the following matrix exactly equivalent to a diagonal matrix  ?

(({{1, 2}, {-2, -3}}))

8.13	There is an example due to Wilkinson, which shows that in the following tridiagonal matrix

(({{10, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {1, 9, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 1, 8, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 1, 7, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 1, 6, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 1, 5, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 1, 4, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 1, 3, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 1, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 1, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 1, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 4, 1, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 5, 1, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 6, 1, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 7, 1, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 8, 1, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 9, 1}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 10}}))

(it is a 21 x 21 matrix, whose main diagonal is 10, 9, ..., 1, 0, 1,...,10, whose parallel diagonals consist of 1's)  the elements in the parallel diagonals need not be small if some of the eigenvalues differ by a small amount only. Prepare this matrix by as few commands as possible and compute the eigenvalues.