E7441: Scientific computing in biology and
biomedicine
Non-square linear systems
Vlad Popovici, Ph.D.
RECETOX
Outline
1 Non-square systems
The underdetermined case
The overdetermined case
2 Numerical methods for LS problem
Orthogonal transformations
Singular Value Decomposition
Total least squares
3 Comparison of various decompositions
4 Eigenvalue problems
Eigenvalue problems
Special forms
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 2 / 57
The systems of linear equations
General form:
Ax = b


a11 . . . a1n
...
am1 . . . amn




x1
...
xn


=


b1
...
bm


if m < n: underdetermined case; find a minimum-norm solution
if m > n: overdetermined case; minimize the squared error
if m = n: determined case; already discussed
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 3 / 57
Reminder
two vectors y, z are orthogonal if yT
z = 0
the span of a set of n independent vectors is span({v1, . . . , vn}) =
n
i=1 αivi | αi ∈ R
the row (column) space of a matrix A is the linear subspace
generated (or spanned) by the rows (colums) of A. Its dimension is
equal to rank(A) ≤ min(m, n).
by definition, span(A) is the column space of A and can be written as
C(A) = {v ∈ Rm
: v = Ax, x ∈ Rn
},
so it is the space of transformed vectors by the action of multiplication
by the matrix.
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 4 / 57
Underdetermined case
m < n there are more variables than equations, hence the solution is
not unique
consider the rows to be linearly independent
then, any n-dimensional vector x ∈ Rn
can be decomposed into
x = x+
+ x−
where x+ is in the row space of A and x− is in the null space of A
(orthogonal to the previous space):
x+
= AT
α Ax−
= 0
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 5 / 57
this leads to
A(x+
+ x−
) = AAT
α + Ax−
= AAT
α = b
AAT
is a m × m nonsingular matrix, so AAT
α = b has a unique
solution α0 = (AAT
)−1
b
the corresponding minimal norm solution to original system is
x+
0
= AT
(AAT
)−1
b
note, however, that the orthogonal component x− remains unspecified
the matrix AT
(AAT
)−1
is called the right pseudo-inverse of A (right:
A · AT
(AAT
)−1
= I)
Python: scipy.linalg.pinv() or numpy.linalg.pinv()
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 6 / 57
Example: let A = [1 2] and b = [3] (hence m = 1).
solution space:
x2 = −
1
2
x1 +
3
2
is a solution, for any x1 ∈ R.
x+ = AT
α =
1
2
α (row space)
Ax− = 0 ⇒ [1 2][x−
1
x−
2
]T
= 0.
⇒ x−
2
= −1
2 x−
1
(null space)
The minimal norm solution is the intersection of solution space with the
row space and is the closest vector to the origin, among all vectors in the
solution space:
x+
0
= [0.6 1.2]T
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 7 / 57
Overdetermined case
if the rows of A are independent, there is no perfect solution to the
system (b span(A))
one needs some other criterion to call a solution acceptable
least squares solution x0 minimizes the square Euclidean norm of the
residual vector:
x0 = arg min
x
∥r∥2
2 = arg min
x
∥b − Ax∥2
2
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 8 / 57
Solution to the LS problem
From a linear system problem, we arrived at solving an optimization
problem with objective function
J =
1
2
∥b − Ax∥2
2 =
1
2
(b − Ax)T
(b − Ax)
Set the derivative wrt x to zero:
∂
∂x
J = AT
b − AT
Ax = 0
which leads to normal equations AT
Ax = AT
b, with the solution
x0 = (AT
A)−1
AT
b
A† = (AT
A)−1
AT
is the left pseudo-inverse of A.
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 9 / 57
Solution to the LS problem - geometric interpretation
let y = Ax, where x is the LS solution
the residual r = b − y is orthogonal to span(A),
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 10 / 57
LS data approximation
Model: y = c3x2
+ c2x + c1. Problem: ci =? when (xi, yi) are given.
See Example 1 in Jupyter notebook.
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 11 / 57
Condition number
if rank(A) = n (columns are independent), the condition number is
cond(A) = ∥A∥2∥A†
∥2
by convention, if rank(A) < n, cond(A) = ∞
for non-square matrices, the condition number measures the
closeness to rank deficiency
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 12 / 57
Numerical methods for LS problem
the LS solution can be obtained using the pseudo-inverse
A† = (AT
A)−1
AT
or by solving the normal equations
AT
Ax = AT
b
which is a system of n equations
AT
A is symmetric positive definite, so it admits a Cholesky
decomposition,
AT
A = LLT
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 13 / 57
Issues with normal equations method
floating-point computations in AT
A and AT
b may lead to information
loss
sensitivity of the solution is worsen, since cond(AT
A) = [cond(A)]2
Example:
Let A =


1 1
ϵ 0
0 ϵ


with ϵ ∈ R+ and ϵ <
√
ϵmach. Then, in floating-point
arithmetic, AT
A =
1 + ϵ2
1
1 1 + ϵ2 =
1 1
1 1
which is singular!
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 14 / 57
Augmented systems
idea: find the solution and the residual as a solution of an extended
system, under the orthogonality requirement
the new system is
I A
AT
0
r
x
=
b
0
despite requiring more storage and not being positive definite, it
allows more freedom in choosing pivots for LU decomposition
in some cases it is useful, but not much used in practice
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 15 / 57
Orthogonal transformations
a matrix Q is orthogonal if QT
Q = I
multiplication of a vector by an orthogonal matrix does not change its
Euclidean norm:
∥Qv∥2
2 = (Qv)T
Qv = vT
QT
Qv = vT
v = ∥v∥2
2
so, multiplying the two sides of the system by Q does not change the
solution
again: try to transform the system so it’s easy to solve e.g. triangular
system
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 16 / 57
an upper triangular overdetermined (m > n) LS problem has the form
R
0
x ≈
b1
b2
where R is an n × n upper triangular matrix and b is partitioned
accordingly
the residual becomes
∥r∥2
2 = ∥b1 − Rx∥2
2 + ∥b2∥2
2
to minimize the residual, one has to minimize ∥b1 − Rx∥2
2
(since ∥b2∥2
2
is fixed) and this leads to the system
Rx = b1
which can be solved by back-substitution
the residual becomes ∥r∥2
2
= ∥b2∥2
2
and x is the LS solution
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 17 / 57
QR factorization
problem: find an m × m orthogonal matrix Q such that an m × n matrix
A can be written as
A = Q
R
0
where R is n × n upper triangular
the new problem to solve is
QT
Ax =
R
0
x ≈
b1
b2
= QT
b
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 18 / 57
if Q is partitioned as Q = [Q1Q2] with Q1 having n columns, then
A = Q
R
0
= Q1R
is called reduced QR factorization of A (Python:
scipy.linalg.qr())
columns of Q1 form an orthonormal basis of span(A), and the
columns of Q2 form an orthonormal basis of span(A)⊥
Q1QT
1
is orthogonal projector onto span(A)
the solution to the initial problem is given by the solution to the square
system
QT
1 Ax = QT
1 b
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 19 / 57
QR factorization
In general, for an m × n matrix A, with m > n, the factorization is
A = QR
and
Q is an orthogonal matrix: QT
Q = I ⇔ Q−1
= QT
R is an upper triangular matrix
solving the normal equations (for LS solution) AT
Ax = AT
b comes to
solving
Rx = QT
b
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 20 / 57
Example
See Example 2 in Jupyter notebook.
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 21 / 57
A statistical perspective
Changing a bit the notation, the linear model is
E[y] = Xβ, Cov(y) = σ2
I
It can be shown that the best linear unbiased estimator is
ˆβ = (XT
X)−1
XT
y = R−1
QT
y
for a decomposition X = QR. Then ˆy = QQT
y. (Gauss-Markov thm.: LS
estimator has the lowest variance among all unbiased linear estimators.)
Also,
Var(ˆβ) = (XT
X)−1
σ2
= (RT
R)−1
σ2
where σ2
= ∥y − ˆy∥2
/(m − n − 1).
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 22 / 57
Computing the QR factorization
similarly to LU factorization, we nullify entries under the diagonal,
column by column
now, use orthogonal transformations:
▶ Householder transformations
▶ Givens rotations
▶ Gram-Schmidt orthogonalization
Python:scipy.linalg.qr()
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 23 / 57
Householder transformations
H = I − 2
vvT
vT v
, v 0
H is orthogonal and symmetric: H = HT
= H−1
v are chosen such that for a vector a:
Ha =


α
0
...
0


= α


1
0
...
0


= αe1
this leads to v = a − αe1 with α = ±∥a∥2, where the sign is chosen to
avoid cancellation
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 24 / 57
Householder QR factorization
apply, the Householder transformation to nuliffy the entries below
diagonal
the process is applied to each column (of the n) and produces a
transformation of the form
Hn . . . H1A =
R
0
where R is n × n upper triangular
then take Q = H1 . . . Hn
note that the multiplication of H with a vector u is much cheaper than
a general matrix-vector multiplication:
Hu = I − 2
vvT
vT v
u = u − 2
vT
u
vT v
v
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 25 / 57
Gram-Schmidt orthogonalization
idea: given two vectors a1 and a2, we seek orthonormal vectors q1
and q2 having the same span
method: subtract from a2 its projection on a1 and normalize the
resulting vectors
apply this method to each column of A to obtain the classical
Gram-Schmidt procedure
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 26 / 57
Algorithm: Classical Gram-Schmidt
for k = 1 to n do
qk ← ak ;
for j = 1 to k − 1 do
rjk ← qT
j
ak ;
qk ← qk − rjk qj;
rkk ← ∥qk ∥2;
qk ← qk /rkk ;
The resulting matrices Q (with qk as columns) and R (with elements rjk )
form the reduced QR factorization of A.
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 27 / 57
Further topics on QR factorization
if rank(A) < n then R is singular and there are multiple solutions x;
choose the x with the smallest norm
in limited precision, the rank can be lower than the theoretical one,
leading to highly sensitive solutions → an alternative could be the
SVD method (next)
there exists a version, QR with pivoting, that chooses everytime the
column with largest Euclidean norm for reduction → improves stability
in rank deficient scenarios
another method of factorization: Givens rotations - makes one 0 at a
time
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 28 / 57
Singular Value Decomposition - SVD
SVD of an m × n matrix A has the form
A = UΣVT
where U is m × m orthogonal matrix, V is n × n orthogonal matrix, and
Σ is m × n diagonal matrix, with
σii =



0 if i j
σi ≥ 0 if i = j
σi are usually ordered such that σ1 ≥ · · · ≥ σn and are called singular
values of A
the columns ui and vi are called left and right singular vectors of A,
respectively
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 29 / 57
minimum norm solution to Ax ≈ b is
x =
σi 0
uT
i
b
σi
vi
for ill-conditioned or rank-deficient problems, the sum should be taken
over "large enough" σ’s: σi≥ϵ . . .
Euclidean norm: ∥A∥2 = maxi{σi}
Euclidean condition number: cond(A) = maxi{σi}
mini{σi}
Rank of A : rank(A) = #{σi > 0}
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 30 / 57
Pseudoinverse (again)
the pseudoinverse of an m × n matrix A with SVD decomposition
A = UΣVT
is
A+
= VΣ−1
UT
where
[Σ−1
]ii =



1/σi for σi > 0
0 otherwise
pseudoinverse always exists and minimum norm solution to Ax ≈ b is
x = A+b
if A is square and nonsingular, A−1
= A+
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 31 / 57
SVD and subspaces relevant to A
ui for which σi > 0 form the orthonormal basis of span(A)
ui for which σi = 0 form the orthonormal basis of the orthogonal
complement of span(A)
vi for which σi = 0 form the orthonormal basis of the null space of A
vi for which σi > 0 form the orthonormal basis of the orthogonal
complement of the null space of A
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 32 / 57
SVD and matrix approximation
A can be re-written as
A = UΣVT
= σ1u1vT
1 + · · · + σnunvT
n
let Ei = uivT
i
; Ei has rank 1 and requires only m + n storage locations
Eix multiplication requires only m + n multiplications
assuming σ1 ≥ σ2 ≥ . . . σn then by using the largest k singular
values, one obtains the closes approximation of A of rank k:
A ≈
k
i=1
σiEi
many applications to image processing, data compression,
cryptography, etc.
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 33 / 57
Example - image compression
Python: scipy.linalg.svd()
Original image and its approximations using 1,2,3,4,5 and 10 terms:
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 34 / 57
Example - image compression
Python: scipy.linalg.svd()
Original image and its approximations using 1,2,3,4,5 and 10 terms:
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 34 / 57
Example - image compression
Python: scipy.linalg.svd()
Original image and its approximations using 1,2,3,4,5 and 10 terms:
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 34 / 57
Example - image compression
Python: scipy.linalg.svd()
Original image and its approximations using 1,2,3,4,5 and 10 terms:
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 34 / 57
Example - image compression
Python: scipy.linalg.svd()
Original image and its approximations using 1,2,3,4,5 and 10 terms:
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 34 / 57
Example - image compression
Python: scipy.linalg.svd()
Original image and its approximations using 1,2,3,4,5 and 10 terms:
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 34 / 57
Example - image compression
Python: scipy.linalg.svd()
Original image and its approximations using 1,2,3,4,5 and 10 terms:
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 34 / 57
Total least squares
Ax ≊ b
ordinary least squares applies when the error affects only b
what if there is error (uncertainty) in A as well?
total least squares minimizes the orthogonal distances, rather than
vertical distances, between model and data
can be computed using SVD of [A, b]
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 35 / 57
Comparison: work effort
computing AT
A requires about n2
m/2 multiplications and solving the
resulting symmetric system, about n3
/6 multiplications
LS problem solution by Householder QR requires about mn2
− n3
/3
multiplications
if m ≫ n, Householder method requires about twice as much work
normal eqs.
cost of SVD is ≈ (4 . . . 10) × (mn2
+ n3
) depending on implementation
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 36 / 57
Comparison: precision
relative error for normal eqs. is ∼ [cond(A)]2
; if cond(A) ≈ 1/
√
ϵmach,
Cholesky factorization will break
Householder method has a relative error
∼ cond(A) + ∥r∥2[cond(A)]2
which is the best achievable for LS problems
Householder method breaks (in back-substitution step) for
cond(A) ⪅ 1/ϵmach
while Householder method is more general and more accurate than
normal equations, it may not always be worth the additional cost
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 37 / 57
Comparison: precision, cont’d
for (nearly) rank-deficient problems, the pivoting Householder method
produces useful solution, while normal equations method fails
SVD is more precise and more robust than Householder method, but
much more expensive computationally
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 38 / 57
Eigenvalue problems
Standard eigenvalue problem
Given a square matrix A ∈ Mn×n(R), find a scalar λ and a vector
x ∈ Rn
, x 0, such that
Ax = λx.
λ is called eigenvalue and x is called eigenvector
a similar "left" eigenvector can be defined as yT
A = λyT
, but this
would be equivalent to a "right" eigenvalue problem (as above) with
AT
as matrix
the definition can be extended to complex-valued matrices
λ can be complex, even if A ∈ Mn×n(R)
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 39 / 57
Characteristic polynomial
previous eq. is equivalent to (A − λI)x = 0 which admits nonzero
solutions if and only if (A − λI) is singular, i.e.
det(A − λI) = 0
det(. . . ) is the characteristic polynomial of matrix A and its roots λi
are the eigenvalues of A
(from Fundamental Theorem of Algebra) for an n × n matrix there are
n eigenvalues (may not all be real or distinct)
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 40 / 57
reciprocal: a polynomial p(λ) = c0 + c1λ + cn−1λn−1
+ λn
has a
companion matrix


0 0 . . . 0 −c0
1 0 . . . 0 −c1
...
...
...
...
...
0 0 . . . 1 −cn−1


the characteristic polynomial is not used in numerical computation,
because:
▶ finding its roots may imply an infinite number of steps
▶ of the sensitivity of the coefficients
▶ too much work to compute the coefficients and find the roots
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 41 / 57
Example
Let A =
0 1
0 −1
. The characteristic equation is
det(A − λI) = 0 ⇔
λ2
+ λ = 0
with solutions λ1 = 0 and λ2 = −1. For eigenvectors v1, v2 (non-null!):
(A − λ1I)v1 =
0 1
0 −1
v11
v21
=
v21
−v21
:=
0
0
so v21 = 0. We choose v11 such that ∥v1∥ = 1, so v11 = 1. Similarly, for
λ2 = −1 we get v2 =
1/
√
2
−1/
√
2
.
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 42 / 57
See Example 4.
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 43 / 57
Sensitivity of the characteristic polynomial
let A =
1 ϵ
ϵ 1
with ϵ > 0 and slightly smaller than ϵmach
the exact eigenvalues are 1 + ϵ and 1 − ϵ
in floating-point arithmetic,
det(A − λI) = λ2
− 2λ + (1 − ϵ2
) = λ2
− 2λ + 1
with the solution 1 (double root)
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 44 / 57
a simple eigenvalue is a simple solution of the characteristic
polynomial (multiplicity of the root is 1)
a defective matrix has eigenvalues with multiplicity larger than 1,
meaning less than n independent eigenvectors
a nondefective matrix has exactly n linearly independent eigenvectors
and can be diagonalized
Q−1
AQ = Λ
where Q is a nonsingular matrix of eigenvectors
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 45 / 57
Eigen-decomposition
it follows that if A admits n independent eigenvectors, it can be
decomposed (factorized) as
A = QΛQ−1
with Q having the eigenvectors of A as columns, and Λ a diagonal
matrix with eigenvalues on the diagonal
theoretically, A−1
= QΛ−1
Q−1
(if λi 0 and all eigenvalues are
distinct)
if A is normal (AH
A = AH
A) then Q becomes unitary
if A is real symmetric, then Q is orthogonal
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 46 / 57
Eigenvectors
the eigenvectors can be arbitrarily scaled
usually, the eigenvectors are normalized, ∥x∥ = 1
the eigenspace is Sλ = {x|Ax = λx}
a subspace S ⊂ Rn
is invariant if AS ⊆ S
for xi eigenvectors, span({xi}) is an invariant subspace
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 47 / 57
Some useful properties
det(A) = N
i=1 λni
i
, where ni is the multiplicity of eigenvalue λi
tr(A) = N
i=1 niλi
the eigenvalues of A−1
are λ−1
i
(for λi 0)
the eigenvectors of A−1
are the same as those of A
A admits an eigen-decomposition if all eigenvalues are distinct
if A is invertible it does not imply that it can be eigen-decomposed;
reciprocally, if A admits an eigen-decomposition, it does not imply it
can be inverted
A can be inverted if and only if λi 0, ∀i
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 48 / 57
Before solving an eigenvalue problem...
do I need all the eigenvalues?
do I need the eigenvectors as well?
is A real or complex?
is A small, dense or large and sparse?
is there anything special about A? e.g.: symmetric, diagonal,
orthogonal, Hermitian, etc etc
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 49 / 57
Conditioning of EV problems
conditioning of EV problem is different than conditioning of linear
systems for the same matrix
sensitivity is "not uniform" among eigenvectors/eigenvalues
for a simple eigenvalue λ, the condition is 1/∥yH
x|, where x and y are
the corresponding right and left normalized eigenvectors (and yH
is
the conjugate transpose)
so the condition is 1/ cos(x, y)
a perturbation of order ϵ in A may perturb the eigenvalue λ by as
much as ϵ/ cos(x, y)
for special cases of A, special forms of conditioning can be derived
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 50 / 57
Computation - general ideas
a matrix B is similar to A if there exists a nonsingular matrix T such
that B = T−1
AT
if y is an eigenvector of B then x = Ty is an eigenvector of A and
HOMEWORK: prove that A and B have the same eigenvalues
transformations:
▶ shift: A ← A − σI
▶ inversion: A ← A−1
(if A is nonsingular)
▶ power: A ← Ak
▶ polynomial: let p be a polynomial, then A ← p(A)
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 51 / 57
Forms attainable by similarity
For a matrix A with given property, the matrices T and B exist such that
B = T−1
AT has the desired property:
A T B
distinct eigenvalues nonsingular diagonal
real symmetric orthogonal real diagonal
complex Hermitian unitary real diagonal
normal unitary diagonal
arbitrary real orthogonal real block triangular (Schur)
arbitrary unitary upper triangular (Schur)
arbitrary nonsingular almost diagonal
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 52 / 57
If A is diagonal...
the eigenvalues are the diagonal entries
the eigenvectors are the columns of the identity matrix
If a matrix is not diagonalizable, one can obtain a Jordan form:


λ1 1
λ1 1
λ1
λ2 1
λ2
λ3
...
λk 1
λk


Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 53 / 57
If A is triangular (Schur form, in general)...
eigenvalues are the elements on the diagonal
eigenvectors are obtained as follows:
If
A − λI =


U11 u U13
0 0 vT
0 0 U33


is triangular, then U11y = u can be solved for y, so that
x =


y
−1
0


is the corresponding eigenvector
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 54 / 57
Symmetric matrices - Jacobi method
idea: start with a symetric matrix A0 and iteratively form
Ak+1 = JT
k
Ak Jk , where Jk is a plane rotation chosen to annihilate a
symmetric pair of entries in Ak with the goal of diagonalizing A
a rotation matrix has the form
cos θ sin θ
− sin θ cos θ
the problem is to find θ
for A =
a b
b c
and requiring that JT
AJ is diagonal, we obtain
1 + tan θ
a − c
b
− tan2
θ = 0
from which we use the root with the smallest magnitude
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 55 / 57
for more general matrices, there are other methods like Power
iterations, with or without deflation, etc.
a generalized eigenvalue problem,
Ax = λBx,
can be solved using the QZ algorithm
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 56 / 57
Singular Value Decomposition - again
we saw that SVD of a m × n matrix A has the form
A = UΣVT
where U is m × m orthogonal matrix and V is n × n orthogonal matrix
and Σ is m × n diagonal matrix with non-negative elements on the
diagonal
this is a eigenvalue-like problem
the columns of U and V are the left and right singular vectors,
respectively and σii are the singular values
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 57 / 57
The relation between SVD and the eigen-decomposition
SVD can be applied to any m × n matrix, while the
eigen-decomposition is applied only to square matrices
the singular values are non-negative while the eigenvalues can be
negative
let A = UΣVT
be SVD of A ⇒ AT
A = (VΣT
UT
)(UΣVT
) = VΣT
ΣVT
also, AT
A is symmetric real matrix, so it has a eigendecomposition
AT
A = QΛQT
, with Q orthogonal. By unicity of decompositions, it
follows that
ΣT
Σ = Λ
V = Q
so σi =
√
λi
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 58 / 57
Questions?
Vlad Popovici, Ph.D. (RECETOX) E7441: Scientific computing in biology and biomedicine 59 / 57