%-----------------------------------------------------
%from Higham&Higham, MATLAB Guide, Ch.9 Linear Algebra (SIAM)
%-----------------------------------------------------
% (1) Norms and Condition Numbers
x=1:4;
[norm(x,1) norm(x,2) norm(x,inf) norm(x,-inf)]%:=min_i|x_i|
A=[1 2 3;4 5 6;7 8 9];
[norm(A,1) norm(A,2) norm(A,inf) norm(A,'fro')]
spy(A)
A=gallery('grcar',8);
[cond(A,1) 1/rcond(A) condest(A)]
%run twise; condest result may vary!
%-----------------------------------------------------
% (2) Linear Equations A*X=B => X=A\B
%uses of \
% ----------square A:
% general A : LU, partial pivoting ("warning off")
% if A = T : solved by substitution
% if A = A*, pos.def.: Cholesky factorization use
X = hilb(15)\ones(15,1);
% ----------mXn A; m>n, overdetermined:
% least squares solution: min norm(A*x-b)
% if r<n, gives basic solution (free var.=0)
B = [1 -1 1;2 3 4; 3 2 5; 4 1 6]; % m=4, n=3,; r=2
b=ones(3,1)
x=B\b;
x_min=pinv(B)*b % pseudoinverse; gives min_norm solution
A = gallery('lauchli',3,0.25), b=[1 2 4 8]';
x = A\b; % least squares solution
xn = lsqnonneg(A,b); % Nonnegative least squares solution
[x xn], [norm(A*x-b) norm(A*xn-b)]
% ----------mXn A; m<n, underdetermined:
C=[1 1 1;1 1 -1], b = [3;1]
x=C\b; % basic solution (free = 0)
y=pinv(A)*b; % minimal solution
[x y]
[norm(x) norm(y)]
%-----------------------------------------------------
% (3) Inverse, Pseudoinverse and Determinant
A = pascal(3), X = inv(A),
norm(A*X - eye(3))
pinv(ones(3))
Y = pinv(A)
norm(A*Y-eye(3))
B = vander(1:5), det(B)
%-----------------------------------------------------
% (4) LU and Cholesky factorization
format short g
A = gallery('fielder',3), [L,U] = lu(A), b=ones(3,1);
x=A\b
y=U\(L\b)
B = pascal(4), R=chol(B), R'*R
[R,p]=chol(B)
%-----------------------------------------------------
% (4) QR factorization
format short e, A
[Q,R] = qr(A)
A(2,2) = eps;
[Q,R,P] = qr(A); R, P
% help qrdelete, qrinsert, qrupdate
%-----------------------------------------------------
% (5) SVD
A = reshape(1:9,3,3);format short e
svd(A)
format
rank(A)
null(A)
orth(A)
rref(A) %recuced row echelon form--very sensitive to roundoff
% A = UCX*, B = VSX*, C* C + S* S = I
[U,V,X,C,S]=gsvd(A,B)
%-----------------------------------------------------
% (6) Eigenvalues and Generalized Eigenvalues
% eig(A), eig(A,B)
[V,D] = eig([2 -1;1 0])
[V,D] = eig([2 -1;-1 1]) % if hermitian, evals sorted, V unitary
norm(V*V'-eye(2))
F = gallery('frank',5)
e=eig(F)
1./e
[V,D,s]=condeig(A)
[diag(D)';s']
A=gallery('triw',3), B=magic(3)
[V,D] = eig(A,B); V, eivals=diag(D)'
A=gallery('fiedler',3), B=gallery('moller',3)
format short g
[V,D] = eig(A,B); V, eivals = diag(D)'
V'*A*V
V'*B*V
% Gen Schur decomp QAZ = T, QBZ = S
% (Q,Z unitary; T,S upper-triang)
% [T,S,Q,Z,V] = qz(A,B)
A=eye(2);B=[20 -10;-10 20];C=[15 -5;-5 15];
[X,e] = polyeig(C,B,A) ;(r^2A + rB +C)x = 0
%-----------------------------------------------------
% (8) Iterations
%---Linear equations
A=gallery('wathen',12,12);n=length(A)
b=ones(n,1);
x=pcg(A,b);
x=pcg(A,b,1e-6,100);
[x,flag,relres,iter] = pcg(A,b,1e-6,100,diag(diag(A)));
% bicg Biconjugate gradient
% bicgstab " " stabilized
% cgs Conjugate " squared
% lsqr " " on normal equations
% pcg (Hpd) prec. " "
% gmres Gener. minimum residual
% minres (H) " "
% qmr quasi " "
% symmlq (H) Symmetric LQ
%---Eigensolvers
A = delsq(numgrid('N',40));
n = length(A);
mnz(A)/n^2; % Percentage of nonzeros
tic, e_all = eig(full(A))';toc
e_all(n:-1:n-4)
options.disp=0; % Turn off intermediate output
tic,e_big = eigs(A,5,'LA',options)';
%-----------------------------------------------------
% (9) Matrix functions
% funm(A,'function') (expm, sqrtm, logm)
A=[17 8 1 0;8 18 8 1;1 8 18 8 0 1 8 17]
sqrtm(A)
A = gallery('jordblock',4,1)
X = expm(A)
real(logm(X))
%-----------------------------------------------------
% (10) Sparse matrices
A = sparse([1 2 2 4 4],[3 1 4 2 4],1:5)
B = full(A)
nnz(A)
whos
A = sparse(i,j,s,m,n)
A = sparse(i,j,s,m,n,nzmax)
[i,j,s]=find(A);
[m,n] = size(A);
B=[1 2 0;1 2 3;0 2 3;0 2 3]
d=[-2 0 1]
A=spdiags(B,d,4,4)
full(A)
n=5;e=ones(n,1);
A=spdiags([-e 4*e -e],[-1 0 1],n,n);
full(A)
% A=sprand(m,n,density,rc)
% lu(A,thresh) : 0, none; 1,partial
% (in between: piv>=thresh*max)
%A=gallery('wathen',8,8);
A=gallery('grcar',8);
subplot(1,2,1), spy(A), subplot(1,2,2), spy(chol(A))
figure
r=symrcm(A);
subplot(1,2,1), spy(A(r,r)), subplot(1,2,2), spy(chol(A(r,r)))
figure
m=symmmd(A);
subplot(1,2,1), spy(A(m,m)), subplot(1,2,2), spy(chol(A(m,m)))
% Viewing the structure of a matrix "spy"
%-----------------------------------------------------