Previous Next Contents

9  Linear Algebra

9.1   aff2ab linear (affine) function to A,b conversion




CALLING SEQUENCE :

[A,b]=aff2ab(afunction,dimX,D [,flag])



PARAMETERS :




DESCRIPTION :

aff2ab returns the matrix representation of an affine function (in the canonical basis).

afunction is a function with imposed syntax: Y=afunction(X,D) where X=list(X1,X2,...,Xp) is a list of p real matrices, and Y=list(Y1,...,Yq) is a list of q real real matrices which depend linearly of the Xi's. The (optional) input D contains parameters needed to compute Y as a function of X. (It is generally a list of matrices).

dimX is a p x 2 matrix: dimX(i)=[nri,nci] is the actual number of rows and columns of matrix Xi. These dimensions determine na, the column dimension of the resulting matrix A: na=nr1*nc1 +...+ nrp*ncp.

If the optional parameter flag='sp' the resulting A matrix is returned as a sparse matrix.

This function is useful to solve a system of linear equations where the unknown variables are matrices.




EXAMPLE :

// Lyapunov equation solver (one unknown variable, one constraint)
deff('Y=lyapunov(X,D)','[A,Q]=D(:);Xm=X(:); Y=list(A''*Xm+Xm*A-Q)')
A=rand(3,3);Q=rand(3,3);Q=Q+Q';D=list(A,Q);dimX=[3,3];
[Aly,bly]=aff2ab(lyapunov,dimX,D);
[Xl,kerA]=linsolve(Aly,bly); Xv=vec2list(Xl,dimX); lyapunov(Xv,D)
Xm=Xv(:); A'*Xm+Xm*A-Q

// Lyapunov equation solver with redundant constraint X=X'
// (one variable, two constraints) D is global variable
deff('Y=ly2(X,D)','[A,Q]=D(:);Xm=X(:); Y=list(A''*Xm+Xm*A-Q,Xm''-Xm)')
A=rand(3,3);Q=rand(3,3);Q=Q+Q';D=list(A,Q);dimX=[3,3];
[Aly,bly]=aff2ab(ly2,dimX,D);
[Xl,kerA]=linsolve(Aly,bly); Xv=vec2list(Xl,dimX); ly2(Xv,D)

// Francis equations
// Find matrices X1 and X2 such that:
// A1*X1 - X1*A2 + B*X2 -A3 = 0
// D1*X1 -D2 = 0 
deff('Y=bruce(X,D)','[A1,A2,A3,B,D1,D2]=D(:),...
[X1,X2]=X(:);Y=list(A1*X1-X1*A2+B*X2-A3,D1*X1-D2)')
A1=[-4,10;-1,2];A3=[1;2];B=[0;1];A2=1;D1=[0,1];D2=1;
D=list(A1,A2,A3,B,D1,D2);
[n1,m1]=size(A1);[n2,m2]=size(A2);[n3,m3]=size(B);
dimX=[[m1,n2];[m3,m2]];
[Af,bf]=aff2ab(bruce,dimX,D);
[Xf,KerAf]=linsolve(Af,bf);Xsol=vec2list(Xf,dimX)
bruce(Xsol,D)

// Find all X which commute with A
deff('y=f(X,D)','y=list(D(:)*X(:)-X(:)*D(:))')
A=rand(3,3);dimX=[3,3];[Af,bf]=aff2ab(f,dimX,list(A));
[Xf,KerAf]=linsolve(Af,bf);[p,q]=size(KerAf);
Xsol=vec2list(Xf+KerAf*rand(q,1),dimX);
C=Xsol(:); A*C-C*A



See Also : linsolve X

9.2   balanc matrix or pencil balancing




CALLING SEQUENCE :

[Ab,X]=balanc(A)
[Eb,Ab,X,Y]=balanc(E,A)



PARAMETERS :




DESCRIPTION :

Balance a square matrix to improve its condition number.

[Ab,X] = balanc(A) finds a similarity transformation X such that Ab = inv(X)*A*X has approximately equal row and column norms.

For matrix pencils,balancing is done for improving the generalized eigenvalue problem.

[Eb,Ab,X,Y] = balanc(E,A) returns left and right transformations X and Y such that Eb=X*E*Y Ab=X*A*Y


REMARK :

Balancing is made in the functions bdiag and spec.


EXAMPLE :

A=[1/2^10,1/2^10;2^10,2^10];
[Ab,X]=balanc(A);
norm(A(1,:))/norm(A(2,:))
norm(Ab(1,:))/norm(Ab(2,:))



See Also : bdiag X

9.3   bdiag block diagonalization, generalized eigenvectors




CALLING SEQUENCE :

[Ab [,X [,bs]]]=bdiag(A [,rmax])



PARAMETERS :




DESCRIPTION :

[Ab [,X [,bs]]]=bdiag(A [,rmax])
performs the block-diagonalization of matrix A. bs gives the structure of the blocks (respective sizes of the blocks). X is the change of basis i.e Ab = inv(X)*A*X is block diagonal.

rmax controls the conditioning of X; the default value is the l1 norm of A.

To get a diagonal form (if it exists) choose a large value for rmax (rmax=1/%eps for example). Generically (for real random A) the blocks are (1x1) and (2x2) and X is the matrix of eigenvectors.


EXAMPLE :

//Real case: 1x1 and 2x2 blocks
a=rand(5,5);[ab,x,bs]=bdiag(a);ab
//Complex case: complex 1x1 blocks
[ab,x,bs]=bdiag(a+%i*0);ab



See Also : schur X, sylv X, spec X

9.4   chfact sparse Cholesky factorization




CALLING SEQUENCE :

spcho=chfact(A)



PARAMETERS :




DESCRIPTION :

spcho=chfact(A) computes the sparse Cholesky factors of sparse matrix A, assumed symmetric positive definite. This function is based on the Ng-Peyton programs (ORNL). See the Fortran programs for a complete description of the variables in spcho. This function is to be used with chsolve.


See Also : chsolve X, sparse X, lufact X, luget X, spchol X

9.5   chol Cholesky factorization




CALLING SEQUENCE :

[R]=chol(X)



PARAMETERS  :




DESCRIPTION :

If X is positive definite, then R = chol(X) produces an upper triangular matrix R such that R'*R = X.

chol(X) uses only the diagonal and upper triangle of X. The lower triangular is assumed to be the (complex conjugate) transpose of the upper.


EXAMPLE :

W=rand(5,5)+%i*rand(5,5);
X=W*W';
R=chol(X);
norm(R'*R-X)



See Also : spchol X, qr X, svd X, bdiag X, fullrf X

9.6   chsolve sparse Cholesky solver




CALLING SEQUENCE :

sol=chsolve(spcho,rhs)



PARAMETERS :




DESCRIPTION :

sol=chsolve(spcho,rhs) computes the solution of sol=A*rhs, with A a symmetric sparse positive definite matrix. This function is based on the Ng-Peyton programs (ORNL). See the Fortran programs for a complete description of the variables in spcho.


EXAMPLE :

A=sprand(20,20,0.1);
A=A*A'+eye();  
spcho=chfact(A);
sol=(1:20)';rhs=A*sol;
spcho=chfact(A);
chsolve(spcho,rhs)



See Also : chfact X, sparse X, lufact X, luget X, spchol X

9.7   coff resolvent (cofactor method)




CALLING SEQUENCE :

[N,d]=coff(M [,var])



PARAMETERS :




DESCRIPTION :

coff computes R=(s*eye()-M)^-1 for M a real matrix. R is given by N/d.

N = numerator polynomial matrix.

d = common denominator.

var character string ('s' if omitted)


EXAMPLE :

M=[1,2;0,3];
[N,d]=coff(M)
N/d
inv(%s*eye()-M)



See Also : coffg X, ss2tf X, nlev X, poly X

9.8   colcomp column compression, kernel, nullspace




CALLING SEQUENCE :

[W,rk]=colcomp(A [,flag] [,tol])



PARAMETERS :




DESCRIPTION :

Column compression of A: Ac = A*W is column compressed i.e

Ac=[0,Af] with Af full column rank, rank(Af) = rank(A) = rk.

flag and tol are optional parameters: flag = 'qr' or 'svd' (default is 'svd').

tol = tolerance parameter (of order %eps as default value).

The ma-rk first columns of W span the kernel of A when size(A)=(na,ma)


EXAMPLE :

A=rand(5,2)*rand(2,5);
[X,r]=colcomp(A);
norm(A*X(:,1:$-r),1)



See Also : rowcomp X, fullrf X, fullrfk X, kernel X


Author : F.D.

9.9   companion companion matrix




CALLING SEQUENCE :

A=companion(p)



PARAMETERS :




DESCRIPTION :

Returns a matrix A with characteristic polynomial equal to p if p is monic. If p is not monic the characteristic polynomial of A is equal to p/c where c is the coefficient of largest degree in p.

If p is a vector of monic polynomials, A is block diagonal, and the characteristic polynomial of the ith block is p(i).


EXAMPLE :

s=poly(0,'s');
p=poly([1,2,3,4,1],'s','c')
det(s*eye()-companion(p))
roots(p)
spec(companion(p))



See Also : spec X, poly X, randpencil X


Author : F.D.

9.10   cond condition number




CALLING SEQUENCE :

cond(X)



PARAMETERS :




DESCRIPTION :

Condition number in 2-norm. cond(X) is the ratio of the largest singular value of X to the smallest.


EXAMPLE :

A=testmatrix('hilb',6);
cond(A)



See Also : rcond X, svd X

9.11   det determinant




CALLING SEQUENCE :

det(X)
[e,m]=det(X)



PARAMETERS :




DESCRIPTION :

det(X) ( m*10^e is the determinant of the square matrix X.

For polynomial matrix det(X) is equivalent to determ(X).

For rational matrices det(X) is equivalent to detr(X).


EXAMPLE :

x=poly(0,'x');
det([x,1+x;2-x,x^2])
w=ssrand(2,2,4);roots(det(systmat(w))),trzeros(w)   //zeros of linear system
A=rand(3,3);
det(A), prod(spec(A))



See Also : detr X, determ X

9.12   ereduc computes matrix column echelon form by qz transformations




CALLING SEQUENCE :

[E,Q,Z [,stair [,rk]]]=ereduc(X,tol)



PARAMETERS :




DESCRIPTION :

Given an m x n matrix X (not necessarily regular) the function ereduc computes a unitary transformed matrix E=Q*X*Z which is in column echelon form (trapezoidal form). Furthermore the rank of matrix X is determined.


EXAMPLE :

X=[1 2 3;4 5 6]
[E,Q,Z ,stair ,rk]=ereduc(X,1.d-15)



See Also : fstair X


Author : Th.G.J. Beelen (Philips Glass Eindhoven). SLICOT

9.13   exp element-wise exponential




CALLING SEQUENCE :

exp(X)



PARAMETERS :




DESCRIPTION :

exp(X) is the (element-wise) exponential of the entries of X.




EXAMPLE :

x=[1,2,3+%i];
log(exp(x))  //element-wise
2^x
exp(x*log(2))



See Also : coff X, log X, expm X

9.14   expm square matrix exponential




CALLING SEQUENCE :

expm(X)



PARAMETERS :




DESCRIPTION :

X is a square matrix expm(X) is the matrix expm(X) = I + X + X^2 /2 + ... The computation is performed by first block-diagonalizing X and then applying a Pade approximation on each block.


EXAMPLE :

X=[1 2;3 4]
expm(X)
logm(expm(X))



See Also : logm X, bdiag X, coff X, log X, exp X

9.15   fstair computes pencil column echelon form by qz transformations




CALLING SEQUENCE :

[AE,EE,QE,ZE,blcks,muk,nuk,muk0,nuk0,mnei]=fstair(A,E,Q,Z,stair,rk,tol)



PARAMETERS :




DESCRIPTION :

Given a pencil sE-A where matrix E is in column echelon form the function fstair computes according to the wishes of the user a unitary transformed pencil QE(sEE-AE)ZE which is more or less similar to the generalized Schur form of the pencil sE-A. The function yields also part of the Kronecker structure of the given pencil.

Q,Z are the unitary matrices used to compute the pencil where E is in column echelon form (see ereduc)


Author : Th.G.J. Beelen (Philips Glass Eindhoven). SLICOT


See Also : quaskro X, ereduc X

9.16   fullrf full rank factorization




CALLING SEQUENCE :

[Q,M,rk]=fullrf(A,[tol])



PARAMETERS :




DESCRIPTION :

Full rank factorization : fullrf returns Q and M such that A = Q*M with range(Q)=range(A) and ker(M)=ker(A), Q full column rank , M full row rank, rk = rank(A) = #columns(Q) = #rows(M).

tol is an optional real parameter (default value is sqrt(%eps)). The rank rk of A is defined as the number of singular values larger than norm(A)*tol.

If A is symmetric, fullrf returns M=Q'.


EXAMPLE :

A=rand(5,2)*rand(2,5);
[Q,M]=fullrf(A);
norm(Q*M-A,1)
[X,d]=rowcomp(A);Y=X';
svd([A,Y(:,1:d),Q])        //span(Q) = span(A) = span(Y(:,1:2))



See Also : svd X, qr X, fullrfk X, rowcomp X, colcomp X


Author : F.D.

9.17   fullrfk full rank factorization of A^k




CALLING SEQUENCE :

[Bk,Ck]=fullrfk(A,k)



PARAMETERS :




DESCRIPTION :

This function computes the full rank factorization of A^k i.e. Bk*Ck=A^k where Bk is full column rank and Ck full row rank. One has range(Bk)=range(A^k) and ker(Ck)=ker(A^k).

For k=1, fullrfk is equivalent to fullrf.


EXAMPLE :

A=rand(5,2)*rand(2,5);[Bk,Ck]=fullrfk(A,3);
norm(Bk*Ck-A^3,1)



See Also : fullrf X, range X


Author : F.D (1990)

9.18   givens Givens transformation




CALLING SEQUENCE :

U=givens(xy)
U=givens(x,y)
[U,c]=givens(xy)
[U,c]=givens(x,y)



PARAMETERS :




DESCRIPTION :

U = givens(x, y) or U = givens(xy) with xy = [x;y] returns a 2x2 unitary matrix U such that:

U*xy=[r;0]=c.

Note that givens(x,y) and givens([x;y]) are equivalent.


EXAMPLE :

A=[3,4;5,6];
U=givens(A(:,1));
U*A



See Also : qr X

9.19   glever inverse of matrix pencil




CALLING SEQUENCE :

[Bfs,Bis,chis]=glever(E,A [,s])



PARAMETERS :




DESCRIPTION :

Computation of (s*E-A)^-1 by generalized Leverrier's algorithm for a matrix pencil. 
(s*E-A)^-1 = (Bfs/chis) - Bis.
chis = characteristic polynomial (up to a multiplicative constant).

Bfs = numerator polynomial matrix.

Bis = polynomial matrix ( - expansion of (s*E-A)^-1 at infinity).

Note the - sign before Bis.




CAUTION  :

This function uses cleanp to simplify Bfs,Bis and chis.


EXAMPLE :

s=%s;F=[-1,s,0,0;0,-1,0,0;0,0,s-2,0;0,0,0,s-1];
[Bfs,Bis,chis]=glever(F)
inv(F)-((Bfs/chis) - Bis)



Author : F. D. (1988)


See Also : rowshuff X, det X, invr X, coffg X, pencan X, penlaur X

9.20   gschur generalized Schur form (matrix pencils).




CALLING SEQUENCE :

[As,Es]=gschur(A,E)
[As,Es,Q,Z]=gschur(A,E)
[As,Es,Z,dim] = gschur(A,E,flag) 
[As,Es,Z,dim]= gschur(A,E,extern)



PARAMETERS :




DESCRIPTION :



Schur form of matrix pencils (QZ algorithm): 
[As,Es] = gschur(A,E)
produces a quasi triangular As matrix and a triangular Es matrix which are the generalized Schur form of the pair A, E. 

[As,Es,Q,Z] = gschur(A,E)
returns in addition two unitary matrices Q and Z such that As=Q*A*Z and Es=Q*E*Z.

Ordered stable form: 
[As,Es,Z,dim] = gschur(A,E,'c')
returns the real generalized Schur form of the pencil s*E-A. In addition, the dim first columns of Z span a basis of the right eigenspace associated with eigenvalues with negative real parts (stable "continuous time" generalized eigenspace). 
[As,Es,Z,dim] = gschur(A,E,'d')
returns the real generalized Schur form of the pencil s*E-A. In addition, the dim first columns of Z make a basis of the right eigenspace associated with eigenvalues with magnitude lower than 1 (stable "discrete time" generalized eigenspace).

General subspace: 
[As,Es,Z,dim] = gschur(A,E,extern)
returns the real generalized Schur form of the pencil s*E-A. In addition, the dim first columns of Z make a basis of the right eigenspace associated with eigenvalues of the pencil which are selected according to a rule which is given by the scilab function extern. (See schur for definition of this function).


EXAMPLE :

s=%s;
F=[-1,s,0,0;0,-1,0,0;0,0,2+s,0;0,0,0,-2+s];
roots(det(F))
[E,A]=pen2ea(F);
[As,Es,Z,dim] = gschur(A,E,'c')
// Other example
a=rand(4,4);b=rand(4,4);[as,bs,qs,zs]=gschur(a,b);
norm(qs*a*zs-as)
norm(qs*b*zs-bs )
clear a;
a(8,8)=2;a(1,8)=1;a(2,[2,3,4,5])=[0.3,0.2,4,6];a(3,[2,3])=[-0.2,.3];
a(3,7)=.5;
a(4,4)=.5;a(4,6)=2;a(5,5)=1;a(6,6)=4;a(6,7)=2.5;a(7,6)=-10;a(7,7)=4;
b=eye(8,8);b(5,5)=0;
[al,be]=gspec(a,b);
[bs,as,q,n]=gschur(b,a,'disc');n-4



See Also : external X, gspec X, pencan X, penlaur X, coffg X, kroneck X

9.21   gspec eigenvalues of matrix pencil




CALLING SEQUENCE :

[al,be]=gspec(A,E)
[al,be,Z]=gspec(A,E)



PARAMETERS :




DESCRIPTION :

[al,be] = gspec(A,E) returns the spectrum of the matrix pencil s E - A, i.e. the roots of the polynomial matrix s E - A. The eigenvalues are given by al./be and if be(i) = 0 the ith eigenvalue is at infinity. (For E = eye(A), al./be is spec(A)).

[al,be,Z] = gspec(A,E) returns in addition the matrix Z of generalized right eigenvectors of the pencil.


EXAMPLE :

A=rand(3,3);
[al,be,Z] = gspec(A,eye(A));al./be
clean(inv(Z)*A*Z)  //displaying the eigenvalues (generic matrix)
A=A+%i*rand(A);E=rand(A);
roots(det(%s*E-A))   //complex case



See Also : gschur X, balanc X, spec X, kroneck X

9.22   hess Hessenberg form




CALLING SEQUENCE :

H = hess(A)
[U,H] = hess(A)



PARAMETERS :




DESCRIPTION :

[U,H] = hess(A) produces a unitary matrix U and a Hessenberg matrix H so that A = U*H*U' and U'*U = Identity. By itself, hess(A) returns H.

The Hessenberg form of a matrix is zero below the first subdiagonal. If the matrix is symmetric or Hermitian, the form is tridiagonal.


EXAMPLE :

A=rand(3,3);[U,H]=hess(A);
and( abs(U*H*U'-A)<1.d-10 )



See Also : qr X, contr X, schur X

9.23   householder Householder orthogonal reflexion matrix




CALLING SEQUENCE :

u=householder(v [,w])



PARAMETERS :




DESCRIPTION :

given 2 column vectors v, w of same size, householder(v,w) returns a unitary column vector u, such that (eye()-2*u*u')*v is proportional to w. (eye()-2*u*u') is the orthogonal Householder reflexion matrix .

w default value is eye(v). In this case vector (eye()-2*u*u')*v is the vector eye(v)*norm(v).


See Also : qr X, givens X

9.24   im_inv inverse image




CALLING SEQUENCE :

[X,dim]=im_inv(A,B [,tol])
[X,dim,Y]=im_inv(A,B, [,tol])



PARAMETERS :




DESCRIPTION :

[X,dim]=im_inv(A,B) computes (A^-1)(B) i.e vectors whose image through A are in range(B)

The dim first columns of X span (A^-1) (B)

tol is a threshold used to test if subspace inclusion; default value is tol = 100*%eps. If Y is returned, then [Y*A*X,Y*B] is partitioned as follows: 
[A11,A12;0,A22],[B1;0]
where B1 has full row rank (equals rank(B)) and A22 has full column rank and has dim columns.


EXAMPLE :

A=[rand(2,5);[zeros(3,4),rand(3,1)]];B=[[1,1;1,1];zeros(3,2)];
W=rand(5,5);A=W*A;B=W*B;
[X,dim]=im_inv(A,B)
svd([A*X(:,1:dim),B])   //vectors A*X(:,1:dim) belong to range(B)
[X,dim,Y]=im_inv(A,B);[Y*A*X,Y*B]



See Also : rowcomp X, spaninter X, spanplus X, linsolve X


Author : F. D.

9.25   inv matrix inverse




CALLING SEQUENCE :

inv(X)



PARAMETERS :




DESCRIPTION :

inv(X) is the inverse of the square matrix X. A warning message is printed if X is badly scaled or nearly singular.

For polynomial matrices or rational matrices in transfer representation, inv(X) is equivalent to invr(X).

For linear systems in state-space representation (syslin list), invr(X) is equivalent to invsyslin(X).


EXAMPLE :

A=rand(3,3);inv(A)*A
//
x=poly(0,'x');
A=[x,1,x;x^2,2,1+x;1,2,3];inv(A)*A
//
A=[1/x,2;2+x,2/(1+x)]
inv(A)*A
//
A=ssrand(2,2,3);
W=inv(A)*A
clean(ss2tf(W))



See Also : slash X, backslash X, pinv X, qr X, lufact X, lusolve X, invr X, coff X, coffg X

9.26   kernel kernel, nullspace




CALLING SEQUENCE :

W=kernel(A [,tol,[,flag])



PARAMETERS :




DESCRIPTION :

W=kernel(A) returns the kernel (nullspace) of A.

flag and tol are optional parameters: flag = 'qr' or 'svd' (default is 'svd').

tol = tolerance parameter (of order %eps as default value).


EXAMPLE :

A=rand(3,1)*rand(1,3);
A*kernel(A)
A=sparse(A);
clean(A*kernel(A))



See Also : colcomp X, fullrf X, fullrfk X, linsolve X


Author : F.D.

9.27   kroneck Kronecker form of matrix pencil




CALLING SEQUENCE :

[Q,Z,Qd,Zd,numbeps,numbeta]=kroneck(F)
[Q,Z,Qd,Zd,numbeps,numbeta]=kroneck(E,A)



PARAMETERS :




DESCRIPTION :

Kronecker form of matrix pencil: kroneck computes two orthogonal matrices Q, Z which put the pencil F=s*E -A into upper-triangular form: 

           | sE(eps)-A(eps) |        X       |      X     |      X        |
           |----------------|----------------|------------|---------------|
           |        O       | sE(inf)-A(inf) |      X     |      X        |
Q(sE-A)Z = |---------------------------------|----------------------------|
           |                |                |            |               |
           |        0       |       0        | sE(f)-A(f) |      X        |
           |--------------------------------------------------------------|
           |                |                |            |               |
           |        0       |       0        |      0     | sE(eta)-A(eta)|
The dimensions of the four blocks are given by:

eps=Qd(1) x Zd(1), inf=Qd(2) x Zd(2), f = Qd(3) x Zd(3), eta=Qd(4)xZd(4)

The inf block contains the infinite modes of the pencil.

The f block contains the finite modes of the pencil

The structure of epsilon and eta blocks are given by:

numbeps(1) = # of eps blocks of size 0 x 1

numbeps(2) = # of eps blocks of size 1 x 2

numbeps(3) = # of eps blocks of size 2 x 3 etc...

numbeta(1) = # of eta blocks of size 1 x 0

numbeta(2) = # of eta blocks of size 2 x 1

numbeta(3) = # of eta blocks of size 3 x 2 etc...

The code is taken from T. Beelen (Slicot-WGS group).


EXAMPLE :

F=randpencil([1,1,2],[2,3],[-1,3,1],[0,3]);
Q=rand(17,17);Z=rand(18,18);F=Q*F*Z;
//random pencil with eps1=1,eps2=1,eps3=1; 2 J-blocks @ infty 
//with dimensions 2 and 3
//3 finite eigenvalues at -1,3,1 and eta1=0,eta2=3
[Q,Z,Qd,Zd,numbeps,numbeta]=kroneck(F);
[Qd(1),Zd(1)]    //eps. part is sum(epsi) x (sum(epsi) + number of epsi) 
[Qd(2),Zd(2)]    //infinity part
[Qd(3),Zd(3)]    //finite part
[Qd(4),Zd(4)]    //eta part is (sum(etai) + number(eta1)) x sum(etai)
numbeps
numbeta



See Also : gschur X, gspec X, systmat X, pencan X, randpencil X, trzeros X

9.28   linsolve linear equation solver




CALLING SEQUENCE :

[x0,kerA]=linsolve(A,b [,x0])



PARAMETERS :




DESCRIPTION :

linsolve computes all the solutions to A*x+b=0.

x0 is a particular solution (if any) and kerA= nullspace of A. Any x=x0+kerA*w with arbitrary w satisfies A*x+b=0.

If compatible x0 is given on entry, x0 is returned. If not a compatible x0, if any, is returned.


EXAMPLE :

A=rand(5,3)*rand(3,8);
b=A*ones(8,1);[x,kerA]=linsolve(A,b);A*x+b   //compatible b
b=ones(5,1);[x,kerA]=linsolve(A,b);A*x+b   //uncompatible b
A=rand(5,5);[x,kerA]=linsolve(A,b), -inv(A)*b  //x is unique



See Also : inv X, pinv X, colcomp X, im_inv X

9.29   lu LU factors of Gaussian elimination




CALLING SEQUENCE :

[L,U]= lu(A)
[L,U,E]= lu(A)



PARAMETERS :




DESCRIPTION :

[L,U]= lu(A) produces two matrices L and U such that A = L*U with U upper triangular and E*L lower triangular for a permutation matrix E.

If A has rank k, rows k+1 to n of U are zero.

 [L,U,E]= lu(A) produces three matrices L, U and E such that E*A = L*U with U upper triangular and E*L lower triangular for a permutation matrix E.




REMARK :

If A is a real matrix, using the function lufact and luget it is possible to obtain the permutation matrices and also when A is not full rank the column compression of the matrix L.

 
[h,rk]=lufact(sparse(a))  // lufact works with sparse real matrices 
[P,L,U,Q]=luget(h)
ludel(h)
P=full(P);L=full(L);U=full(U);Q=full(Q); 
// P,Q are permutation matrices P*L*U*Q=A



See Also : lufact X, luget X, lusolve X, qr X, svd X

9.30   ludel utility function used with lufact




CALLING SEQUENCE :

ludel(hand)



PARAMETERS :




DESCRIPTION :

This function is used in conjunction with lufact. It clears the internal memory space used to store the result of lufact.

The sequence of commands [p,r]=lufact(A);x=lusolve(p,b);ludel(p); solves the sparse linear system A*x = b and clears p.


See Also : sparse X, lufact X, luget X

9.31   lufact sparse lu factorization




CALLING SEQUENCE :

[hand,rk]=lufact(A,prec)



PARAMETERS :




DESCRIPTION :

[hand,rk]=lufact(A) performs the lu factorization of sparse matrix A. hand (no display) is used by lusolve (for solving linear system) and luget (for retrieving the factors). hand should be cleared by the command: ludel(hand);

The A matrix needs not be full rank but must be square (since A is assumed sparse one may add zeros if necessary to squaring down A).


EXAMPLE :

a=rand(5,5);b=rand(5,1);A=sparse(a);
[h,rk]=lufact(A);
x=lusolve(h,b);a*x-b
ludel(h)



See Also : sparse X, lusolve X, luget X

9.32   luget sparse lu factorization




CALLING SEQUENCE :

[P,L,U,Q]=luget(ptr)



PARAMETERS :




DESCRIPTION :

[P,L,U,Q]=luget(ptr) with ptr obtained by the command [ptr,rk]=lufact(A) with A a sparse matrix returns four sparse matrices such that P*L*U*Q=A.

The A matrix needs not be full rank but must be square (since A is assumed sparse one may add zeros if necessary to squaring down A).

If A is singular, the L matrix is column compressed (with rk independent nonzero columns): the nonsingular sparse matrix Q'*inv(U) column compresses A.


EXAMPLE :

a=rand(5,2)*rand(2,5);A=sparse(a);
[ptr,rk]=lufact(A);[P,L,U,Q]=luget(ptr);
full(L), P*L*U*Q-A
clean(P*L*U*Q-A)
ludel(ptr)



See Also : sparse X, lusolve X, luget X, clean X

9.33   lusolve sparse linear system solver




CALLING SEQUENCE :

lusolve(hand,b)
lusolve(A,b)



PARAMETERS :




DESCRIPTION :

x=lusolve(hand,b) solves the sparse linear system A*x = b.

[hand,rk]=lufact(A) is the output of lufact.

x=lusolve(A,b) solves the sparse linear system \fVA*x = b\fR.


EXAMPLE :

non_zeros=[1,2,3,4];rows_cols=[1,1;2,2;3,3;4,4];
sp=sparse(rows_cols,non_zeros);
[h,rk]=lufact(sp);x=lusolve(h,[1;1;1;1]);ludel(h)
rk,sp*x


non_zeros=[1,2,3,4];rows_cols=[1,1;2,2;3,3;4,4];
sp=sparse(rows_cols,non_zeros);
x=lusolve(sp,-ones(4,1));
sp*x



See Also : sparse X, lufact X, slash X, backslash X

9.34   lyap Lyapunov equation




CALLING SEQUENCE :

[X]=lyap(A,C,'c')
[X]=lyap(A,C,'d')



PARAMETERS :




DESCRIPTION :

X = lyap(A,C,flag) solves the continuous time or discrete time matrix Lyapunov matrix equation: 
 
A'*X + X*A = C       ( flag='c' )
A'*X*A - X = C       ( flag='d' )
Note that a unique solution exist if and only if an eigenvalue of A is
not an eigenvalue of -A (flag='c')  or 1 over an eigenvalue of A 
(flag='d').



EXAMPLE :

A=rand(4,4);C=rand(A);C=C+C';
X=lyap(A,C,'c');
A'*X + X*A -C
X=lyap(A,C,'d');
A'*X*A - X -C



See Also : sylv X, ctr_gram X, obs_gram X

9.35   nlev Leverrier's algorithm




CALLING SEQUENCE :

[num,den]=nlev(A,z [,rmax])



PARAMETERS :




DESCRIPTION :

[num,den]=nlev(A,z [,rmax]) computes: 
(z*eye()-A)^(-1)
by block diagonalization of A followed by Leverrier's algorithm on each block.


REMARK :

This algorithm is better than the usual leverrier algorithm but still not perfect!


EXAMPLE :

A=rand(3,3);x=poly(0,'x');
[NUM,den]=nlev(A,'x')
clean(den-poly(A,'x'))
clean(NUM/den-inv(x*eye()-A))



See Also : coff X, coffg X, glever X, ss2tf X


Author : F. D., S. S.

9.36   orth orthogonal basis




CALLING SEQUENCE :

Q=orth(A)



PARAMETERS :




DESCRIPTION :

Q=orth(A) returns Q, an orthogonal basis for the span of A. Range(Q) = Range(A) and Q'*Q=eye.

The number of columns of Q is the rank of A as determined by the QR algorithm.


EXAMPLE :

A=rand(5,3)*rand(3,4);
[X,dim]=rowcomp(A);X=X';
svd([orth(A),X(:,1:dim)])



See Also : qr X, rowcomp X, colcomp X, range X

9.37   pbig eigen-projection




CALLING SEQUENCE :

[Q,M]=pbig(A,thres,flag)



PARAMETERS :




DESCRIPTION :

Projection on eigen-subspace associated with eigenvalues with real part >= thres (flag='c') or with magnitude >= thres (flag='d').

The projection is defined by Q*M, Q is full column rank, M is full row rank and M*Q=eye.

If flag='c', the eigenvalues of M*A*Q = eigenvalues of A with real part >= thres.

If flag='d', the eigenvalues of M*A*Q = eigenvalues of A with magnitude >= thres.

If flag='c' and if [Q1,M1] = full rank factorization (fullrf) of eye()-Q*M then eigenvalues of M1*A*Q1 = eigenvalues of A with real part < thres.

If flag='d' and if [Q1,M1] = full rank factorization (fullrf) of eye()-Q*M then eigenvalues of M1*A*Q1 = eigenvalues of A with magnitude < thres.


EXAMPLE :

A=diag([1,2,3]);X=rand(A);A=inv(X)*A*X;
[Q,M]=pbig(A,1.5,'d');
spec(M*A*Q)
[Q1,M1]=fullrf(eye()-Q*M);
spec(M1*A*Q1)



See Also : psmall X, projspec X, fullrf X


Author : F. D. (1988)

9.38   pencan canonical form of matrix pencil




CALLING SEQUENCE :

[Q,M,i1]=pencan(Fs)
[Q,M,i1]=pencan(E,A)



PARAMETERS :




DESCRIPTION :

Given the regular pencil Fs=s*E-A, pencan returns matrices Q and M such than M*(s*E-A)*Q is in "canonical" form.

M*E*Q is a block matrix 
 
[I,0;
 0,N]
with N nilpotent and i1 = size of the I matrix above.

M*A*Q is a block matrix: 
[Ar,0;
 0,I]



EXAMPLE :

F=randpencil([],[1,2],[1,2,3],[]);
F=rand(6,6)*F*rand(6,6);
[Q,M,i1]=pencan(F);
W=clean(M*F*Q)
roots(det(W(1:i1,1:i1)))
det(W($-2:$,$-2:$))



See Also : glever X, penlaur X, rowshuff X


Author : F. D.

9.39   penlaur Laurent coefficients of matrix pencil




CALLING SEQUENCE :

[Si,Pi,Di,order]=penlaur(Fs)
[Si,Pi,Di,order]=penlaur(E,A)



PARAMETERS :




DESCRIPTION :

penlaur computes the first Laurent coefficients of (s*E-A)^-1 at infinity.

(s*E-A)^-1 = ... + Si/s - Pi - s*Di + ... at s = infinity.

order = order of the singularity (order=index-1).

The matrix pencil Fs=s*E-A should be invertible.

For a index-zero pencil, Pi, Di,... are zero and Si=inv(E).

For a index-one pencil (order=0),Di =0.

For higher-index pencils, the terms -s^2 Di(2), -s^3 Di(3),... are given by:

Di(2)=Di*A*Di, Di(3)=Di*A*Di*A*Di (up to Di(order)).


REMARK :

Experimental version: troubles when bad conditioning of so*E-A


EXAMPLE :

F=randpencil([],[1,2],[1,2,3],[]);
F=rand(6,6)*F*rand(6,6);[E,A]=pen2ea(F);
[Si,Pi,Di]=penlaur(F);
[Bfs,Bis,chis]=glever(F);
norm(coeff(Bis,1)-Di,1)



See Also : glever X, pencan X, rowshuff X


Author : F. D. (1988,1990)



9.40   pinv pseudoinverse




CALLING SEQUENCE :

pinv(A,[tol])



PARAMETERS :




DESCRIPTION :

X = pinv(A) produces a matrix X of the same dimensions as A' such that:

A*X*A = A, X*A*X = X and both A*X and X*A are Hermitian .

The computation is based on SVD and any singular values lower than a tolerance are treated as zero: this tolerance is accessed by X=pinv(A,tol).


EXAMPLE :

A=rand(5,2)*rand(2,4);
norm(A*pinv(A)*A-A,1)



See Also : rank X, svd X, qr X

9.41   polar polar form




CALLING SEQUENCE :

[Ro,Theta]=polar(A)



PARAMETERS :




DESCRIPTION :

[Ro,Theta]=polar(A) returns the polar form of A i.e.:

A=Ro*expm(%i*Theta) Ro symmetric >=0 and Theta hermitian >=0.


EXAMPLE :

A=rand(5,5);
[Ro,Theta]=polar(A);
norm(A-Ro*expm(%i*Theta),1)



See Also : expm X, svd X


Author : F. D.

9.42   proj projection




CALLING SEQUENCE :

P = proj(X1,X2)



PARAMETERS :




DESCRIPTION :

P is the projection on X2 parallel to X1.


EXAMPLE :

X1=rand(5,2);X2=rand(5,3);
P=proj(X1,X2);
norm(P^2-P,1)
trace(P)    // This is dim(X2)
[Q,M]=fullrf(P);
svd([Q,X2])   // span(Q) = span(X2)



See Also : projspec X, orth X, fullrf X


Author : F. D.

9.43   projspec spectral operators




CALLING SEQUENCE :

[S,P,D,i]=projspec(A)



PARAMETERS :




DESCRIPTION :

Spectral characteristics of A at 0.

S = reduced resolvent at 0 (S = -Drazin_inverse(A)).

P = spectral projection at 0.

D = nilpotent operator at 0.

index = index of the 0 eigenvalue.

One has (s*eye()-A)^(-1) = D^(i-1)/s^i +... + D/s^2 + P/s - S - s*S^2 -... around the singularity s=0.


EXAMPLE :

deff('j=jdrn(n)','j=zeros(n,n);for k=1:n-1;j(k,k+1)=1;end')
A=sysdiag(jdrn(3),jdrn(2),rand(2,2));X=rand(7,7);
A=X*A*inv(X);
[S,P,D,index]=projspec(A);
index   //size of J-block
trace(P)  //sum of dimensions of J-blocks
A*S-(eye()-P)
norm(D^index,1)



See Also : coff X


Author : F. D.

9.44   psmall spectral projection




CALLING SEQUENCE :

[Q,M]=psmall(A,thres,flag)



PARAMETERS :




DESCRIPTION :

Projection on eigen-subspace associated with eigenvalues with real part < thres (flag='c') or with modulus < thres (flag='d').

The projection is defined by Q*M, Q is full column rank, M is full row rank and M*Q=eye.

If flag='c', the eigenvalues of M*A*Q = eigenvalues of A with real part < thres.

If flag='d', the eigenvalues of M*A*Q = eigenvalues of A with magnitude < thres.

If flag='c' and if [Q1,M1] = full rank factorization (fullrf) of eye()-Q*M then eigenvalues of M1*A*Q1 = eigenvalues of A with real part >= thres.

If flag='d' and if [Q1,M1] = full rank factorization (fullrf) of eye()-Q*M then eigenvalues of M1*A*Q1 = eigenvalues of A with magnitude >= thres.


EXAMPLE :

A=diag([1,2,3]);X=rand(A);A=inv(X)*A*X;
[Q,M]=psmall(A,2.5,'d');
spec(M*A*Q)
[Q1,M1]=fullrf(eye()-Q*M);
spec(M1*A*Q1)



See Also : pbig X, proj X, projspec X


Author : F. D. (1988)

9.45   qr QR decomposition




CALLING SEQUENCE :

[Q,R]=qr(X)
[Q,R,E]=qr(X)
[Q,R,rk,E]=qr(X [,tol])



PARAMETERS :




DESCRIPTION :

[Q,R] = qr(X) produces an upper triangular matrix R of the same dimension as X and a unitary matrix Q so that X = Q*R.

[Q,R,E] = qr(X) produces a (column) permutation matrix E, an upper triangular R with decreasing diagonal elements and a unitary Q so that X*E = Q*R.

[Q,R,rk,E] = qr(X ,tol) returns rk = rank estimate of X i.e. rk is the number of diagonal elements in R which are larger than tol.

[Q,R,rk,E] = qr(X) returns rk = rank estimate of X i.e. rk is the number of diagonal elements in R which are larger than R(1,1)*%eps*max(size(R).


EXAMPLE :

A=rand(5,2)*rand(2,5);
[Q,R,rk,E] = qr(A,1.d-10);
norm(Q'*A-R)
svd([A,Q(:,1:rk)])    //span(A) =span(Q(:,1:rk))



See Also : rank X, svd X, rowcomp X, colcomp X

9.46   quaskro quasi-Kronecker form




CALLING SEQUENCE :

[Q,Z,Qd,Zd,numbeps,numbeta]=quaskro(F)
[Q,Z,Qd,Zd,numbeps,numbeta]=quaskro(E,A)
[Q,Z,Qd,Zd,numbeps,numbeta]=quaskro(F,tol)
[Q,Z,Qd,Zd,numbeps,numbeta]=quaskro(E,A,tol)



PARAMETERS :




DESCRIPTION :

Quasi-Kronecker form of matrix pencil: quaskro computes two orthogonal matrices Q, Z which put the pencil F=s*E -A into upper-triangular form: 

           | sE(eps)-A(eps) |        X       |      X     |
           |----------------|----------------|------------|
           |        O       | sE(inf)-A(inf) |      X     |
Q(sE-A)Z = |=================================|============|
           |                                 |            |
           |                O                | sE(r)-A(r) |
The dimensions of the blocks are given by:

eps=Qd(1) x Zd(1), inf=Qd(2) x Zd(2), r = Qd(3) x Zd(3)

The inf block contains the infinite modes of the pencil.

The f block contains the finite modes of the pencil

The structure of epsilon blocks are given by:

numbeps(1) = # of eps blocks of size 0 x 1

numbeps(2) = # of eps blocks of size 1 x 2

numbeps(3) = # of eps blocks of size 2 x 3 etc...

The complete (four blocks) Kronecker form is given by the function kroneck which calls quaskro on the (pertransposed) pencil sE(r)-A(r).

The code is taken from T. Beelen


See Also : kroneck X, gschur X, gspec X

9.47   randpencil random pencil




CALLING SEQUENCE :

F=randpencil(eps,infi,fin,eta)



PARAMETERS :




DESCRIPTION :

Utility function. F=randpencil(eps,infi,fin,eta) returns a random pencil F with given Kronecker structure. The structure is given by: eps=[eps1,...,epsk]: structure of epsilon blocks (size eps1x(eps1+1),....) fin=[l1,...,ln] set of finite eigenvalues (assumed real) (possibly []) infi=[k1,...,kp] size of J-blocks at infinity ki>=1 (infi=[] if no J blocks). eta=[eta1,...,etap]: structure ofeta blocks (size eta1+1)xeta1,...)

epsi's should be >=0, etai's should be >=0, infi's should be >=1.

If fin is a (monic) polynomial, the finite block admits the roots of fin as eigenvalues.

If fin is a vector of polynomial, they are the finite elementary divisors of F i.e. the roots of p(i) are finite eigenvalues of F.


EXAMPLE :

F=randpencil([0,1],[2],[-1,0,1],[3]);
[Q,Z,Qd,Zd,numbeps,numbeta]=kroneck(F);
Qd, Zd
s=poly(0,'s');
F=randpencil([],[1,2],s^3-2,[]); //regular pencil
det(F)



See Also : kroneck X, pencan X, penlaur X

9.48   range range (span) of A^k




CALLING SEQUENCE :

[X,dim]=range(A,k)



PARAMETERS :




DESCRIPTION :

Computation of Range A^k ; the first dim columns of X span the range of A^k.


See Also : fullrfk X, rowcomp X


Author : F. D.

9.49   rank rank




CALLING SEQUENCE :

[i]=rank(X)
[i]=rank(X,tol)



PARAMETERS :




DESCRIPTION :

rank(X) is the numerical rank of X i.e. the number of singular values of X that are larger than norm(size(X),'inf') * norm(X) * %eps.

rank(X,tol) is the number of singular values of X that are larger than tol.


REMARK :

Note that the default value of tol is proportional to norm(X). As a consequence

rank([1.d-80,0;0,1.d-80]) is 2 !.


EXAMPLE :

rank([1.d-80,0;0,1.d-80])
rank([1,0;0,1.d-80])



See Also : svd X, qr X, rowcomp X, colcomp X, lu X

9.50   rcond inverse condition number




CALLING SEQUENCE :

rcond(X)



PARAMETERS :




DESCRIPTION :

rcond(X) is an estimate for the reciprocal of the condition of X in the 1-norm.

If X is well conditioned, rcond(X) is close to 1. If not, rcond(X) is close to 0.

[r,z]=rcond(X) sets r to rcond(X) and returns z such that

norm(X*z,1) = r*norm(X,1)*norm(z,1) Thus, if rcond is small z is a vector in the kernel.


EXAMPLE :

A=diag([1:10]);
rcond(A)
A(1,1)=0.000001;
rcond(A)



See Also : svd X, cond X, inv X

9.51   rowcomp row compression, range




CALLING SEQUENCE :

[W,rk]=rowcomp(A [,flag] [,tol])



PARAMETERS :




DESCRIPTION :

Row compression of A. Ac = W*A is a row compressed matrix: i.e. Ac=[Af;0] with Af full row rank.

flag and tol are optional parameters: flag='qr' or 'svd' (default 'svd').

tol is a tolerance parameter (of order sqrt(%eps) as default value).

The rk first columns of W' span the range of A.

The rk first (top) rows of W span the row range of A.


REMARK :

A non zero vector x belongs to range(A) iff W*x is row compressed in accordance with Ac i.e the norm of its last components is small w.r.t its first components.


EXAMPLE :

A=rand(5,2)*rand(2,4);   // 4 col. vectors, 2 independent.
[X,dim]=rowcomp(A);Xp=X';
svd([Xp(:,1:dim),A])     //span(A) = span(Xp(:,1:dim)
x=A*rand(4,1);      //x belongs to span(A)
y=X*x  
norm(y(dim+1:$))/norm(y(1:dim))    // small



See Also : colcomp X, fullrf X, fullrfk X


Author : F. D.

9.52   rowshuff shuffle algorithm




CALLING SEQUENCE :

[Ws,Fs1]=rowshuff(Fs, [alfa])



PARAMETERS :




DESCRIPTION :

Shuffle algorithm: Given the pencil Fs=s*E-A , returns Ws=W(s) (square polynomial matrix) such that:

Fs1 = s*E1-A1 = W(s)*(s*E-A) is a pencil with non singular E1 matrix.

This is possible iff the pencil Fs = s*E-A is regular (i.e. invertible). The degree of Ws is equal to the index of the pencil.

The poles at infinity of Fs are put to alfa and the zeros of Ws are at alfa.

Note that (s*E-A)^-1 = (s*E1-A1)^-1 * W(s) = (W(s)*(s*E-A))^-1 *W(s)


EXAMPLE :

F=randpencil([],[2],[1,2,3],[]);
F=rand(5,5)*F*rand(5,5);   // 5 x 5 regular pencil with 3 evals at 1,2,3
[Ws,F1]=rowshuff(F,-1);
[E1,A1]=pen2ea(F1);
svd(E1)           //E1 non singular
roots(det(Ws))
clean(inv(F)-inv(F1)*Ws,1.d-7)



See Also : pencan X, glever X, penlaur X


Author : F. D.

9.53   rref computes matrix row echelon form by lu transformations




CALLING SEQUENCE :

R=rref(A)



PARAMETERS :




DESCRIPTION :

rref computes the row echelon form of the given matrix by left lu decomposition. If ones need the transformation used just call X=rref([A,eye(m,m)]) the row echelon form R is X(:,1:n) and the left transformation L is given by X(:,n+1:n+m) such as L*A=R


EXAMPLE :

A=[1 2;3 4;5 6];
X=rref([A,eye(3,3)]);
R=X(:,1:2)
L=X(:,3:5);L*A



See Also : lu X, qr X

9.54   schur [ordered] Schur decomposition




CALLING SEQUENCE :

[U,T] = schur(A) 
[U,dim]=schur(A,flag)
[U,dim]=schur(A,myfunction)



PARAMETERS :




DESCRIPTION :

Schur forms, ordered Schur forms


USUAL SCHUR FORM :

[U,T] = schur(A) produces a Schur matrix T and a unitary matrix U so that A = U*T*U' and U'*U = eye(U). By itself, schur(A) returns T. If A is complex, the Complex Schur Form is returned in matrix T. The Complex Schur Form is upper triangular with the eigenvalues of A on the diagonal. If A is real, the Real Schur Form is returned. The Real Schur Form has the real eigenvalues on the diagonal and the complex eigenvalues in 2-by-2 blocks on the diagonal.


ORDERED STABLE FORM :

[T,dim]=schur(A,'c') returns an unitary matrix T which transforms A into schur form. In addition, the dim first columns of T make a basis of the eigenspace of A associated with eigenvalues with negative real parts (stable "continuous time" eigenspace).

[T,dim]=schur(A,'d') returns an unitary matrix T which transforms A into schur form. In addition, the dim first columns of T span a basis of the eigenspace of A associated with eigenvalues with magnitude lower than 1 (stable "discrete time" eigenspace).


GENERAL EIGENSPACE :

[T,dim]=schur(A,a_function) returns an unitary matrix T which transforms A into schur form. In addition, the dim first columns of T span a basis of the eigenspace of A associated with the eigenvalues which are selected by the function a_function.

This function must be of the following type (here a_function is "rule"): 
function [flag]=rule(x)

flag=...
x is a vector with three components which characterizes either a real eigenvalue or a pair of complex conjugate eigenvalues.

If x(1)=1, a real eigenvalue is considered and this eigenvalue is x(2)/x(3).

If x(1)=2, a pair of complex conjugate eigenvalues is considered. The sum of these two eigenvalues (twice the real part) is x(2) and the product (squared magnitude) is x(3).

On return, flag should be 1 if the real eigenvalue is selected or the pair of eigenvalues is selected and 0 otherwise.


EXAMPLE OF FUNCTION :

        function [flag]=disc(x)
        ls =x(1);flag=0;
        select  ls
           case 1 then if abs(x(2)) < ro*abs(x(3)) then flag=1;end
           case 2 then if x(3) < ro*ro then flag=1;end
        end
The function disc selects the eigenvalues with magnitude lower than a given scalar ro. And for ro=1 the calling sequence [T,dim]=schur(A,'d') and [T,dim]=schur(A,disc) are equivalent.

Another useful example is %choose (see function code in SCIDIR/macros/percent)


EXAMPLE :

A=diag([-0.9,-2,2,0.9]);X=rand(A);A=inv(X)*A*X;
[U,d]=schur(A,'c');
A1=U'*A*U;
spec(A1(1:d,1:d))      //stable cont. eigenvalues
[U,d]=schur(A,'c');
A1=U'*A*U;
spec(A1(1:d,1:d))      //stable disc. eigenvalues



See Also : gschur X, ricc X, pbig X, psmall X

9.55   spaninter subspace intersection




CALLING SEQUENCE :

[X,dim]=spaninter(A,B [,tol])



PARAMETERS :




DESCRIPTION :

[X,dim]=spaninter(A,B) computes the intersection of range(A) and range(B).

The first dim columns of X span this intersection i.e. X(:,1:dim) is an orthogonal basis for range(A) inter range(B)

In the X basis A and B are respectively represented by:

X'*A and X'*B.

tol is a threshold (sqrt(%eps) is the default value).


EXAMPLE :

A=rand(5,3)*rand(3,4);     // A is 5 x 4, rank=3
B=[A(:,2),rand(5,1)]*rand(2,2);
[X,dim]=spaninter(A,B);
X1=X(:,1:dim);    //The intersection
svd(A),svd([X1,A])   // X1 in span(A)
svd(B),svd([B,X1])   // X1 in span(B)



See Also : spanplus X, spantwo X


Author : F. D.

9.56   spanplus sum of subspaces




CALLING SEQUENCE :

[X,dim,dima]=spanplus(A,B[,tol])



PARAMETERS :




DESCRIPTION :

[X,dim,dima]=spanplus(A,B) computes a basis X such that:

the first dima columns of X span Range(A) and the following (dim-dima) columns make a basis of A+B relative to A.

The dim first columns of X make a basis for A+B.

One has the following canonical form for [A,B]: 
         [*,*]    (dima rows)
X'*[A,B]=[0,*]    (dim-dima rows)
         [0,0]
tol is an optional argument (see function code).


EXAMPLE :

A=rand(6,2)*rand(2,5);      // rank(A)=2
B=[A(:,1),rand(6,2)]*rand(3,3);   //two additional independent vectors
[X,dim,dimA]=spanplus(A,B);
dimA
dim



See Also : spaninter X, im_inv X, spantwo X


Author : F. D.

9.57   spantwo sum and intersection of subspaces




CALLING SEQUENCE :

[Xp,dima,dimb,dim]=spantwo(A,B, [tol])



PARAMETERS :




DESCRIPTION :

Given two matrices A and B with same number of rows, returns a square matrix Xp (non singular but not necessarily orthogonal) such that : 

         [A1, 0]    (dim-dimb rows)
Xp*[A,B]=[A2,B2]    (dima+dimb-dim rows)
         [0, B3]    (dim-dima rows)
         [0 , 0]
The first dima columns of inv(Xp) span range(A).

Columns dim-dimb+1 to dima of inv(Xp) span the intersection of range(A) and range(B).

The dim first columns of inv(Xp) span range(A)+range(B).

Columns dim-dimb+1 to dim of inv(Xp) span range(B).

Matrix [A1;A2] has full row rank (=rank(A)). Matrix [B2;B3] has full row rank (=rank(B)). Matrix [A2,B2] has full row rank (=rank(A inter B)). Matrix [A1,0;A2,B2;0,B3] has full row rank (=rank(A+B)).


EXAMPLE :

A=[1,0,0,4;
   5,6,7,8;
   0,0,11,12;
   0,0,0,16];
B=[1,2,0,0]';C=[4,0,0,1]; 
Sl=ss2ss(syslin('c',A,B,C),rand(A));
[no,X]=contr(Sl('A'),Sl('B'));CO=X(:,1:no);  //Controllable part
[uo,Y]=unobs(Sl('A'),Sl('C'));UO=Y(:,1:uo);  //Unobservable part
[Xp,dimc,dimu,dim]=spantwo(CO,UO);    //Kalman decomposition
Slcan=ss2ss(Sl,inv(Xp));



See Also : spanplus X, spaninter X


Author : F. D.

9.58   spchol sparse cholesky factorization




CALLING SEQUENCE :

[R,P] = spchol(X)



PARAMETERS :




DESCRIPTION :

[R,P] = spchol(X) produces a lower triangular matrix R such that P*R*R'*P' = X.


EXAMPLE :

X=[
3.,  0.,  0.,  2.,  0.,  0.,  2.,  0.,  2.,  0.,  0. ;
0.,  5.,  4.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0. ;
0.,  4.,  5.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0. ;
2.,  0.,  0.,  3.,  0.,  0.,  2.,  0.,  2.,  0.,  0. ;
0.,  0.,  0.,  0. , 5.,  0.,  0.,  0.,  0.,  0.,  4. ;
0.,  0.,  0.,  0.,  0.,  4.,  0.,  3.,  0.,  3.,  0. ;
2.,  0.,  0.,  2.,  0.,  0.,  3.,  0.,  2.,  0.,  0. ;
0.,  0.,  0.,  0.,  0.,  3.,  0.,  4.,  0.,  3.,  0. ;
2.,  0.,  0.,  2.,  0.,  0.,  2.,  0.,  3.,  0.,  0. ;
0.,  0.,  0.,  0.,  0.,  3.,  0.,  3.,  0.,  4.,  0. ;
0.,  0.,  0.,  0.,  4.,  0.,  0.,  0.,  0.,  0.,  5.];
X=sparse(X);[R,P] = spchol(X);
max(P*R*R'*P'-X)



See Also : sparse X, lusolve X, luget X, chol X

9.59   spec eigenvalues




CALLING SEQUENCE :

evals=spec(A)



PARAMETERS :




DESCRIPTION :

evals=spec(A) returns in vector evals the eigenvalues of A.

Eigenvectors are obtained by bdiag.


EXAMPLE :

A=diag([1,2,3]);X=rand(3,3);A=inv(X)*A*X;
spec(A)
//
x=poly(0,'x');
pol=det(x*eye()-A)
roots(pol)
//
[Ab,X,bs]=bdiag(A);
Ab
clean(inv(X)*A*X)



See Also : poly X, det X, gspec X, schur X, bdiag X, colcomp X

9.60   sqroot W*W' hermitian factorization




CALLING SEQUENCE :

sqroot(X)



PARAMETERS :




DESCRIPTION :

W=sqroot(X) returns W such that X=W*W' (uses SVD).


EXAMPLE :

X=rand(5,2)*rand(2,5);X=X*X';
W=sqroot(X)
norm(W*W'-X,1)
//
X=rand(5,2)+%i*rand(5,2);X=X*X';
W=sqroot(X)
norm(W*W'-X,1)



See Also : chol X, svd X

9.61   sva singular value approximation




CALLING SEQUENCE :

[U,s,V]=sva(A,k)
[U,s,V]=sva(A,tol)



PARAMETERS :




DESCRIPTION :

Singular value approximation.

[U,S,V]=sva(A,k) with k an integer >=1, returns U,S and V such that B=U*S*V' is the best L2 approximation of A with rank(B)=k.

[U,S,V]=sva(A,tol) with tol a real number, returns U,S and V such that B=U*S*V' such that L2-norm of A-B is at most tol.


EXAMPLE :

A=rand(5,4)*rand(4,5);
[U,s,V]=sva(A,2);
B=U*s*V';
svd(A)
svd(B)
clean(svd(A-B))



See Also : svd X

9.62   svd singular value decomposition




CALLING SEQUENCE :

s=svd(X)
[U,S,V]=svd(X)
[U,S,V]=svd(X,0)
[U,S,V,rk]=svd(X [,tol])



PARAMETERS :




DESCRIPTION :

[U,S,V] = svd(X) produces a diagonal matrix S , of the same dimension as X and with nonnegative diagonal elements in decreasing order, and unitary matrices U and V so that X = U*S*V'.

[U,S,V] = svd(X,0) produces the "economy size" decomposition. If X is m-by-n with m > n, then only the first n columns of U are computed and S is n-by-n.

s = svd(X) by itself, returns a vector s containing the singular values.

[U,S,V,rk]=svd(X,tol) gives in addition rk, the numerical rank of X i.e. the number of singular values larger than tol.

The default value of tol is the same as in rank.


EXAMPLE :

X=rand(4,2)*rand(2,4)
svd(X)
sqrt(spec(X*X'))



See Also : rank X, qr X, colcomp X, rowcomp X, sva X, spec X

9.63   sylv Sylvester equation.




CALLING SEQUENCE :

sylv(A,B,C,flag)



PARAMETERS :




DESCRIPTION :

X = sylv(A,B,C,'c') computes X, solution of the "continuous time" Sylvester equation 
A*X+X*B=C
X=sylv(A,B,C,'d') computes X, solution of the "discrete time" Sylvester equation 
A*X*B-X=C



EXAMPLE :

A=rand(4,4);C=rand(4,3);B=rand(3,3);
X = sylv(A,B,C,'c');
norm(A*X+X*B-C)
X=sylv(A,B,C,'d') 
norm(A*X*B-X-C)



See Also : lyap X

9.64   trace trace




CALLING SEQUENCE :

trace(X)



PARAMETERS :




DESCRIPTION :

trace(X) is the trace of the matrix X.

Same as sum(diag(X)).


EXAMPLE :

A=rand(3,3);
trace(A)-sum(spec(A))



See Also : det X

Previous Next Contents