Lecture 9: Matrix decompositions and how we compute them

Syllabus

Week 1: Matrices, vectors, matrix/vector norms, scalar products & unitary matrices
Week 2: TAs-week (Strassen, FFT, a bit of SVD)
Week 3: Matrix ranks, singular value decomposition, linear systems, eigenvalues
Week 4: Matrix decompositions: QR, LU, SVD + test + structured matrices start

Recap of the previous lecture

  • Eigenvectors and eigenvalues
  • Characteristic polynomial and why it is a bad idea
  • Power method: $x := Ax$
  • Gershgorin theorem
  • Schur theorem: $A = U T U^*$
  • Normal matrices: $A^* A = A A^*$
  • A brief whiteboard illustration of how to compute the Schur form: QR-algorithm

Today lecture

Today we will talk about matrix factorizations

  • LU decomposition and Gaussian elimination
  • QR decomposition and Gram-Schmidt algorithm
  • Schur decomposition and QR-algorithm

We already introduced it some time ago, but now we are going to do it in more details.

These are the basic matrix factorizations in numerical linear algebra.

General concept of matrix factorization

In numerical linear algebra we need to solve different tasks, for example:

  • Solve linear systems $Ax = f$
  • Compute eigenvalues / eigenvectors
  • Compute singular values / singular vectors
  • Compute inverses, even sometimes determinants
  • Compute matrix functions like $\exp(A), \cos(A)$.

In order to do this, we replace the matrix as a sum and/or product of matrices with simpler structure,
such that we can solve abovementioned tasks faster / more stable form.

What is a simpler structure?

What is a simpler structure

We already encountered several classes of matrices with structure.

In dense matrix business, the most important classes are

unitary matrices ($U^* U = I$, why?), upper/lower triangular matrices , diagonal matrices

Other classes of matrices

For sparse matrices the sparse constraints are often included in the factorizations.

For Toeplitz matrices an important class of matrices is the matrices with low displacement rank,

which is based on the low-rank matrices.

The class of low-rank matrices and block low-rank matrices is everywhere.

Plan

The plan for todays lecture is to discuss the decompositions one-by-one and point out:

  • How to compute a particular decomposition
  • When the decomposition exists
  • What is done in "real life"

Decompositions we want to discuss today

  • LU factorization
  • Positive definite matrices and Cholesky factorization
  • QR-decomposition and Gram-Schmidt algorithm
  • 1 slide about the SVD (more tomorrow)

What is the LU factorization

LU factorization is the representation of a given matrix $A$ as a product
of a lower triangular matrix $L$ and an upper triangular matrix $U$:

$$A = LU$$

This factorization is non-unique, so it is typical to require that the matrix $L$ has ones on the diagonal.

Does it always exist?

Main goal of the LU decomposition is to solve linear systems, because

$$ A^{-1} f = (L U)^{-1} f = U^{-1} L^{-1} f, $$

and this reduces to the solution of two linear systems $$ L y = f, $$ and $$ U x = y. $$

Solving linear systems with triangular matrices is easy: you just solve the equation with one unknown, put into the second equation, and so on.

The computational cost of such algorithm is $\mathcal{O}(n^2)$ arithmetic operations, once $LU$ factorization is known.

Note that the inverse of triangular matrix is a triangular matrix.

In [2]:
%matplotlib inline
import numpy as np
import sympy
from sympy.matrices import Matrix
import IPython
from sympy.interactive.printing import init_printing
init_printing(use_latex=True)
n = 5
w = Matrix(n, n, lambda i, j: 1/(i + j +  sympy.Integer(1)/2))
L, U, tmp = w.LUdecomposition()
#Generate the final expression
#fn = u'%s \\times %s = %s' % (sympy.latex(L), sympy.latex(U), sympy.latex(w))
#L * U - w
L
⎡ 1    0    0   0   0⎤
⎢                    ⎥
⎢1/3   1    0   0   0⎥
⎢                    ⎥
⎢1/5  6/7   1   0   0⎥
⎢                    ⎥
⎢          15        ⎥
⎢1/7  5/7  ──   1   0⎥
⎢          11        ⎥
⎢                    ⎥
⎢     20   210  28   ⎥
⎢1/9  ──   ───  ──  1⎥
⎣     33   143  15   ⎦

The classical LU-decomposition algorithm

The simplest way to implement LU-decomposition is the following 3-cycle

In [3]:
from sympy.interactive.printing import init_printing
init_printing(use_latex=True)
L = Matrix(n, n, lambda i, j: 0)
U = Matrix(n, n, lambda i, j: 0)
a = w.copy() #Our matrix
for k in xrange(n): #Eliminate one row
    L[k, k] = 1
    for i in xrange(k+1, n):
        L[i, k] = a[i, k] / a[k, k]
        for j in xrange(k+1, n):
            a[i, j] = a[i, j] - L[i, k] * a[k, j]
    for j in xrange(k, n):
        U[k, j] = a[k, j]

Existence of the LU-decomposition

The LU-decomposition algorithm does not fail

if we do not divide by zero at every step.

When it is so, for which class of matrices?

It is true for strictly regular matrices.

Strictly regular matrices and LU decomposition

A matrix $A$ is called strictly regular, if all of its principal minors
(i.e, submatrices in the first $k$ rows and $k$ columns) are non-singular.

In this case, there always exists a LU-decomposition. The reverse is also true (check!).

Proof. If there is a LU-decomposition, then the matrix is strictly regular

This follows from the fact that to get a minor you multiply a corresponding submatrix of $L$ by a corresponding submatrix of $U$, and they are non-singular (since non-singularity of triangular matrices is equivalent to the fact that their diagonal elements are not equal to zero).

The other way can be proven by induction. Suppose that we know that for all $(n-1) \times (n-1)$ matrices will non-singular minors LU-decomposition exists.

Then, consider the block form $$ A = \begin{bmatrix} a & c^{\top} \\ b & D \end{bmatrix}, $$ where $D$ is $(n-1) \times (n-1)$.

Then we do "block elimination":

$$ \begin{bmatrix} 1 & 0 \\ -z & I \end{bmatrix} \begin{bmatrix} a & c^{\top} \\ b & D \end{bmatrix}= \begin{bmatrix} a & c^{\top} \\ 0 & A_1 \end{bmatrix}, $$

where $z = \frac{b}{a}, \quad A_1 = D - \frac{1}{a} b c^{\top}$.
We can show that $A_1$ is also strictly regular, thus it has (by induction) the LU-decomposition:
$$ A_1 = L_1 U_1, $$ And that gives the $LU$ decomposition of the original matrix.

When LU fails

What happens, if the matrix is not strictly regular (or the pivots in the Gaussian elimination are really small?).

There is classical $2 \times 2$ example of a matrix with a bad LU decomposition.

The matrix we look at is
$$ A = \begin{pmatrix} \varepsilon & 1 \\ 1 & 1 \end{pmatrix} $$

If $\varepsilon$ is sufficiently small, we might fail.

Let us do some demo here.

In [45]:
import numpy as np
eps = 1.12e-16
a = [[eps, 1],[1.0,  1]]
a = np.array(a)
a0 = a.copy()
n = a.shape[0]
L = np.zeros((n, n))
U = np.zeros((n, n))
for k in xrange(n): #Eliminate one row   
    L[k, k] = 1
    for i in xrange(k+1, n):
        L[i, k] = a[i, k] / a[k, k]
        for j in xrange(k+1, n):
            a[i, j] = a[i, j] - L[i, k] * a[k, j]
    for j in xrange(k, n):
        U[k, j] = a[k, j]

print 'L * U - A:', np.dot(L, U) - a0
L
L * U - A: [[ 0.  0.]
 [ 0.  0.]]
Out[45]:
array([[  1.00000000e+00,   0.00000000e+00],
       [  8.92857143e+15,   1.00000000e+00]])

The concept of pivoting

We can do pivoting, i.e. permute rows and columns to maximize $A_{kk}$ that we divide over.
The simplest but effective strategy is the row pivoting: at each step, select the index that is maximal in modulus, and put it onto the diagonal.

It gives us the decomposition

$$A = P L U,$$

where $P$ is a permutation matrix.

What makes row pivoting good?

It is made good by the fact that

$$ | L_{ij}|<1, $$

but the elements of $U$ can grow, up to $2^n$! (in practice, this is very rarely encountered).

Can you come up with a matrix where the elements of $U$ grow as much as possible.

Positive definite matrices

Strictly regular matrices have LU-decomposition.

An important subclass of strictly regular matrices is the class of Hermitian positive definite matrices

Positive definite matrices(2)

Definition: A matrix is called positive definite if for any $x: \Vert x \Vert \ne 0$ we have

$$(Ax, x) > 0.$$

Symmetric positive definite matrices

Statement: A Hermitian positive definite matrix $A$ is strictly regular and has Cholesky factorization of the form

$$A = RR^*,$$

where $R$ is upper triangular.

Let us try to prove this fact (on the whiteboard).

The Cholesky factorization is always stable.

It is sometimes referred to as "square root" of the matrix.

Summary on the LU for dense matrices

  • Principal submatrices non-singular -- there exists LU
  • Pivoting is needed to avoid error growth
  • Complexity is $\mathcal{O}(n^3)$ for the factorization step (can be speeded by using Strassen ideas) and $\mathcal{O}(n^2)$ for a solve
  • For SPD matrices Cholesky factorization is the method of choice (stability, twice less memory to store)

QR-decomposition

The next decomposition: QR decomposition. Again from the name it is clear that a matrix is represented as a product
$$ A = Q R, $$ where $Q$ is an orthogonal (unitary) matrix and $R$ is upper triangular.

The matrix sizes: $Q$ is $n \times m$, $R$ is $m \times m$.

QR-factorization is defined for any rectangular matrix.

QR-decomposition: applications

This algorithm plays a crucial role in many problems:

  • Computing orthogonal bases in a linear space
  • Used in the preprocessing step for the SVD
  • QR-algorithm for the computation of eigenvectors and eigenvalues (one of the 10 most important algorithms of the 20th century) is based on the QR-decomposition

Theorem

Every rectangular $n \times m$ matrix, $n \geq m$ matrix has a QR-decomposition. There are several ways to prove it and compute it.

  • (theoretical) Using the Gram matrices and Cholesky factorization
  • (geometrical) Using the Gram-Schmidt orthogonalization
  • (practical) Using Householder/Givens transformations

Mathematical way

If we have the representation of the form $$A = QR,$$ then $A^* A = R^* R$, the matrix $A^* A$ is called Gram matrix, and its elements are scalar products of the columns of $A$.

Math way: full column rank

Assume that $A$ has full column rank. Then,

It is simple to show that $A^* A$ is positive definite:

$$ (A^* A y, y) = (Ay, Ay)^2 = \Vert Ay \Vert \geq 0. $$

Therefore, $A^* A = R^* R$ always exists.

Then the matrix $A R^{-1}$ is unitary:

$$ (A R^{-1})^* (AR^{-1})= R^{-*} A^* A R^{-1} = R^{-*} R^* R R^{-1} = I. $$

Rank-deficient case

When an $n \times m$ matrix does not have full column rank, it is said to be rank-deficient.

The QR-decomposition, however, also exists.

For any rank-deficient matrix there is a sequence of full-column rank matrices $A_k$ such that $A_k \rightarrow A$.

Each $A_k$ can be decomposed as $A_k = Q_k R_k$.

The set of all unitary matrices is compact, thus there exists a converging subsequence $Q_{n_k} \rightarrow Q$,

and $Q^* A_k \rightarrow Q^* A = R$, which is triangular.

QR-decomposition and Gram matrices

So, the simplest way to compute QR-decomposition is then

$$A^* A = R^* R,$$

(Cholesky factorization)

$$Q = A R^{-1}.$$

It is a bad idea for numerical stability. Let us do some demo (for a submatrix of the Hilbert matrix).

In [47]:
import numpy as np
n = 20
r = 8
#a = np.random.randn(n, r)
a = [[1.0/(i+j+0.5) for i in xrange(r)] for j in xrange(n)]
a = np.array(a)
q, Rmat = np.linalg.qr(a)
e = np.eye(r)
print 'Built-in QR orth', np.linalg.norm(np.dot(q.T, q) - e)
gram_matrix = a.T.dot(a)
Rmat1 = np.linalg.cholesky(gram_matrix)
q1 = np.dot(a, np.linalg.inv(Rmat1.T))
print 'Via Gram matrix:', np.linalg.norm(np.dot(q1.T, q1) - e)
Built-in QR orth 8.77419464402e-16
Via Gram matrix: 0.265317882364

Second way: Gram-Schmidt orthogonalization

QR-decomposition is a mathematical way of writing down the Gram-Schmidt orthogonalization process.
Given a sequence of vectors $a_1, \ldots, a_m$ we want to find orthogonal basis $q_1, \ldots, q_m$ such that every $a_i$ is a linear combination of such vectors.

Gram-Schmidt:

  1. $q_1 := a_1/\Vert a_1 \Vert$
  2. $q_2 := a_2 - (a_2, q_1) q_1, \quad q_2 := q_2/\Vert q_2 \Vert$
  3. $q_3 := a_3 - (a_3, q_1) q_1 - (a_3, q_2) q_2, \quad q_2 := q_3/\Vert q_3 \Vert$
  4. And go on

Note that the transformation from $Q$ to $A$ has triangular structure, since from the $k$-th vector we subtract only the previous ones.

Gram-Schmidt is unstable

Gram-Schmidt can be very unstable (i.e., the produced vectors will be not orthogonal, especially if $q_k$ has small norm).
This is called loss of orthogonality.

There is a remedy, called modified Gram-Schmidt method. Instead of doing

$$q_k := a_k - (a_k, q_1) q_1 - \ldots - (a_k, q_{k-1}) q_{k-1}$$

we do it step-by-step. First we set $q_k := a_k$ and orthogonalize sequentially:

$$ q_k := q_k - (q_k, q_1)q_1, \quad q_k := q_{k} - (q_k,q_2)q_2, \ldots $$

In exact arithmetic, it is the same. In floating point it is absolutely different!

Note that the complexity is $\mathcal{O}(nm^2)$ operations

QR-decomposition: the (almost) practical way

If $A = QR$, then
$$ R = Q^* A, $$ and we need to find a certain orthogonal matrix $Q$ that brings a matrix into orthogonal form.
For simplicity, we will look for an $n \times n$ matrix such that $$ Q^* A = \begin{bmatrix}

  • & & \ 0 & & \ 0 & 0 & * \ &0_{(n-m) \times m} \end{bmatrix} $$

We will do it column-by-column.
First, we find a Householder matrix $H_1 = (I - 2 uu^{\top})$ such that (we illustrate on a $4 \times 3$ matrix)

$$ H_1 A = \begin{bmatrix} * & * & * \\ 0 & * & * \\ 0 & * & * \\ 0 & * & * \end{bmatrix} $$

Then, $$ H_2 H_1 A = \begin{bmatrix}

* & * & * \\
0 & * & * \\
0 & 0 & * \\
0 & 0 & *

\end{bmatrix}, $$ where $$ H_2 = \begin{bmatrix} 1 & 0 \\ 0 & H'_2, \end{bmatrix} $$ and $H'_2$ is a $3 \times 3$ Householder matrix.

Finally, $$ H_3 H_2 H_1 A = \begin{bmatrix}

* & * & * \\
0 & * & * \\
0 & 0 & * \\
0 & 0 & 0

\end{bmatrix}, $$

You can try to implement by yourself, it is simple.

QR-decomposition: real life

In reality, since this is a dense matrix factorization, you should implement the algorithm in terms of blocks.

Instead of using Householder transformation, we use block Householder transformation of the form

$$H = (I - 2UU^*), $$

where $U^* U = I$.

This allows us to use BLAS-3 operations.

Rank-revealing QR-decomposition

The QR-decomposition can be also used to compute the (numerical) rank of the matrix.

It is done via so-called rank-revealing factorization.

It is based on the representation

$$P A = Q R, $$

where $P$ is the permutation matrix (it permutes columns), and $R$ has the block form

$$R = \begin{bmatrix} R_{11} & R_{12} \\ 0 & R_{22}\end{bmatrix},$$

and the norm of $R_{22}$ is small, so you can find the numerical rank by looking at it.

Summary

LU and QR decompositions can be computed using direct methods in finite amount of operations.

What about Schur form and SVD?

They can not be computed by direct methods they can only be computed by iterative methods

Schur form

Every matrix can be written as a product

$$A = U T U^*,$$

with triangular $T$ and unitary $U$

and this decomposition gives eigenvalues of the matrix.

QR-algorithm

The QR algorithm (Kublanovskaya and Francis independently proposed it in 1961).

Do not **mix** QR algorithm and QR decomposition.

QR-decomposition is the representation of a matrix, whereas QR algorithm uses QR decomposition computes the eigenvalues!

Way to QR-algorithm

Consider the equation

$$A = Q T Q^*,$$

and rewrite it in the form

$$ Q T = A Q. $$

Derivation of the QR-algorithm as fixed-point

Then we can write down the iterative process

$$ Q_{k+1} R_{k+1} = A Q_k, \quad Q_{k+1}^* A = R_{k+1} Q^*_k $$

Introduce

$$A_k = Q^* _k A Q_k = Q^*_k Q_{k+1} R_{k+1} = \widehat{Q}_k R_{k+1}$$

.

and the new approximation reads

$$A_{k+1} = Q^*_{k+1} A Q_{k+1} = \widehat{Q}^*_k A_k \widehat{Q}_k = R_{k+1} \widehat{Q}_k.$$

Final formulas

The final formulas are then written in the form $$ A_0 = A, \quad A_k = Q_k R_k, \quad A_{k+1} = R_k Q_k. $$

The matrices $A_{k}$ are unitary similar, and $A_k \rightarrow T$.

The complexity of each step is $\mathcal{O}(n^3)$.

In [49]:
import numpy as np
n = 4
a = [[1.0/(i + j + 0.5) for i in xrange(n)] for j in xrange(n)]
niters = 100
for k in xrange(niters):
    q, rmat = np.linalg.qr(a)
    a = rmat.dot(q)
    print 'current a:'
    print a
current a:
[[  2.40047183e+00   1.43485636e-01   4.99605047e-03  -5.56291523e-05]
 [  1.43485636e-01   3.59286592e-01   1.60534687e-02  -1.94225770e-04]
 [  4.99605047e-03   1.60534687e-02   1.60692682e-02  -2.80885670e-04]
 [ -5.56291523e-05  -1.94225770e-04  -2.80885670e-04   2.40684296e-04]]
current a:
[[  2.41031150e+00   2.09437993e-02   3.18722176e-05   5.45279213e-09]
 [  2.09437993e-02   3.50196113e-01   6.87806955e-04   1.28079871e-07]
 [  3.18722176e-05   6.87806955e-04   1.53250849e-02   4.17732654e-06]
 [  5.45279205e-09   1.28079871e-07   4.17732654e-06   2.35678648e-04]]
current a:
[[  2.41051991e+00   3.04114601e-03   2.02621961e-07  -5.33227597e-13]
 [  3.04114601e-03   3.49989111e-01   3.00995528e-05  -8.62074652e-11]
 [  2.02621961e-07   3.00995528e-05   1.53236760e-02  -6.42429491e-08]
 [ -5.33147641e-13  -8.62074827e-11  -6.42429491e-08   2.35677492e-04]]
current a:
[[  2.41052431e+00   4.41545673e-04   1.28806658e-09   1.32059800e-16]
 [  4.41545673e-04   3.49984720e-01   1.31786000e-06   5.80333420e-14]
 [  1.28806664e-09   1.31786000e-06   1.53236733e-02   9.88053919e-10]
 [  5.21260158e-17   5.80510077e-14   9.88053906e-10   2.35677492e-04]]
current a:
[[  2.41052440e+00   6.41081268e-05   8.18816607e-12  -7.99356474e-17]
 [  6.41081268e-05   3.49984627e-01   5.77009674e-08  -2.14109052e-17]
 [  8.18822358e-12   5.77009674e-08   1.53236733e-02  -1.51962427e-11]
 [ -5.09637195e-21  -3.90911829e-17  -1.51962302e-11   2.35677492e-04]]
current a:
[[  2.41052440e+00   9.30787458e-06   5.19949264e-14   7.99300814e-17]
 [  9.30787458e-06   3.49984626e-01   2.52637036e-09  -1.76560777e-17]
 [  5.20524342e-14   2.52637031e-09   1.53236733e-02   2.33729939e-13]
 [  4.98273388e-25   2.63237618e-20   2.33717423e-13   2.35677492e-04]]
current a:
[[  2.41052440e+00   1.35141258e-06   2.73389069e-16  -7.99300126e-17]
 [  1.35141258e-06   3.49984625e-01   1.10614263e-10   1.76826923e-17]
 [  3.30896669e-16   1.10614211e-10   1.53236733e-02  -3.60708065e-15]
 [ -4.87162969e-29  -1.77262591e-23  -3.59456477e-15   2.35677492e-04]]
current a:
[[  2.41052440e+00   1.96211922e-07  -5.54040645e-17   7.99300027e-17]
 [  1.96211922e-07   3.49984625e-01   4.84316731e-12  -1.76827548e-17]
 [  2.10350596e-18   4.84311567e-12   1.53236733e-02   6.78001457e-17]
 [  4.76300288e-33   1.19367537e-26   5.52842647e-17   2.35677492e-04]]
current a:
[[  2.41052440e+00   2.84880569e-08  -5.74941943e-17  -7.99300013e-17]
 [  2.84880568e-08   3.49984625e-01   2.12101875e-13   1.76827614e-17]
 [  1.33719609e-20   2.12050235e-13   1.53236733e-02  -1.33661508e-17]
 [ -4.65679822e-37  -8.03813647e-30  -8.50269816e-19   2.35677492e-04]]
current a:
[[  2.41052440e+00   4.13618795e-09  -5.75074807e-17   7.99300011e-17]
 [  4.13618792e-09   3.49984625e-01   9.33601458e-15  -1.76827623e-17]
 [  8.50053870e-23   9.28437503e-15   1.53236733e-02   1.25289581e-17]
 [  4.55296169e-41   5.41283161e-33   1.30771163e-20   2.35677492e-04]]
current a:
[[  2.41052440e+00   6.00534165e-10  -5.75075651e-17  -7.99300010e-17]
 [  6.00534132e-10   3.49984625e-01   4.58145201e-16   1.76827624e-17]
 [  5.40378175e-25   4.06505655e-16   1.53236733e-02  -1.25160822e-17]
 [ -4.45144049e-45  -3.64496748e-36  -2.01125535e-22   2.35677492e-04]]
current a:
[[  2.41052440e+00   8.71917294e-11  -5.75075656e-17   7.99300010e-17]
 [  8.71916971e-11   3.49984625e-01   6.94379277e-17  -1.76827625e-17]
 [  3.43517725e-27   1.77983814e-17   1.53236733e-02   1.25158841e-17]
 [  4.35218299e-49   2.45449866e-39   3.09330281e-24   2.35677492e-04]]
current a:
[[  2.41052440e+00   1.26594160e-11  -5.75075656e-17  -7.99300010e-17]
 [  1.26593838e-11   3.49984625e-01   5.24188280e-17   1.76827625e-17]
 [  2.18373786e-29   7.79281605e-19   1.53236733e-02  -1.25158811e-17]
 [ -4.25513873e-53  -1.65284428e-42  -4.75748755e-26   2.35677492e-04]]
current a:
[[  2.41052440e+00   1.83805126e-12  -5.75075656e-17   7.99300010e-17]
 [  1.83801902e-12   3.49984625e-01   5.16736663e-17  -1.76827625e-17]
 [  1.38819941e-31   3.41199465e-20   1.53236733e-02   1.25158810e-17]
 [  4.16025833e-57   1.11301516e-45   7.31699713e-28   2.35677492e-04]]
current a:
[[  2.41052440e+00   2.66894672e-13  -5.75075656e-17  -7.99300010e-17]
 [  2.66862429e-13   3.49984625e-01   5.16410403e-17   1.76827625e-17]
 [  8.82476614e-34   1.49390251e-21   1.53236733e-02  -1.25158810e-17]
 [ -4.06749357e-61  -7.49497557e-49  -1.12535128e-29   2.35677492e-04]]
current a:
[[  2.41052440e+00   3.87780641e-14  -5.75075656e-17   7.99300010e-17]
 [  3.87458212e-14   3.49984625e-01   5.16396118e-17  -1.76827625e-17]
 [  5.60989272e-36   6.54087989e-23   1.53236733e-02   1.25158810e-17]
 [  3.97679726e-65   5.04707040e-52   1.73078584e-31   2.35677492e-04]]
current a:
[[  2.41052440e+00   5.65775813e-15  -5.75075656e-17  -7.99300010e-17]
 [  5.62551523e-15   3.49984625e-01   5.16395492e-17   1.76827625e-17]
 [  3.56620174e-38   2.86384884e-24   1.53236733e-02  -1.25158810e-17]
 [ -3.88812328e-69  -3.39866613e-55  -2.66194182e-33   2.35677492e-04]]
current a:
[[  2.41052440e+00   8.49012831e-16  -5.75075656e-17   7.99300010e-17]
 [  8.16769928e-16   3.49984625e-01   5.16395465e-17  -1.76827625e-17]
 [  2.26702996e-40   1.25390319e-25   1.53236733e-02   1.25158810e-17]
 [  3.80142654e-73   2.28864085e-58   4.09405604e-35   2.35677492e-04]]
current a:
[[  2.41052440e+00   1.50829928e-16  -5.75075656e-17  -7.99300010e-17]
 [  1.18587025e-16   3.49984625e-01   5.16395464e-17   1.76827625e-17]
 [  1.44114809e-42   5.49007053e-27   1.53236733e-02  -1.25158810e-17]
 [ -3.71666295e-77  -1.54115666e-61  -6.29664209e-37   2.35677492e-04]]
current a:
[[  2.41052440e+00   4.94605823e-17  -5.75075656e-17   7.99300010e-17]
 [  1.72176791e-17   3.49984625e-01   5.16395464e-17  -1.76827625e-17]
 [  9.16136026e-45   2.40376408e-28   1.53236733e-02   1.25158810e-17]
 [  3.63378941e-81   1.03780540e-64   9.68421076e-39   2.35677492e-04]]
current a:
[[  2.41052440e+00   3.47427422e-17  -5.75075656e-17  -7.99300010e-17]
 [  2.49983902e-18   3.49984625e-01   5.16395464e-17   1.76827625e-17]
 [  5.82386518e-47   1.05246038e-29   1.53236733e-02  -1.25158810e-17]
 [ -3.55276376e-85  -6.98851769e-68  -1.48942780e-40   2.35677492e-04]]
current a:
[[  2.41052440e+00   3.26058554e-17  -5.75075656e-17   7.99300010e-17]
 [  3.62952238e-19   3.49984625e-01   5.16395464e-17  -1.76827625e-17]
 [  3.70222376e-49   4.60807643e-31   1.53236733e-02   1.25158810e-17]
 [  3.47354481e-89   4.70602478e-71   2.29073410e-42   2.35677492e-04]]
current a:
[[  2.41052440e+00   3.22956003e-17  -5.75075656e-17  -7.99300010e-17]
 [  5.26971239e-20   3.49984625e-01   5.16395464e-17   1.76827625e-17]
 [  2.35349899e-51   2.01759313e-32   1.53236733e-02  -1.25158810e-17]
 [ -3.39609227e-93  -3.16900811e-74  -3.52314004e-44   2.35677492e-04]]
current a:
[[  2.41052440e+00   3.22505543e-17  -5.75075656e-17   7.99300010e-17]
 [  7.65110827e-21   3.49984625e-01   5.16395464e-17  -1.76827625e-17]
 [  1.49611634e-53   8.83379887e-34   1.53236733e-02   1.25158810e-17]
 [  3.32036676e-97   2.13399055e-77   5.41857552e-46   2.35677492e-04]]
current a:
[[  2.41052440e+000   3.22440141e-017  -5.75075656e-017  -7.99300010e-017]
 [  1.11086628e-021   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  9.51079275e-056   3.86777697e-035   1.53236733e-002  -1.25158810e-017]
 [ -3.24632976e-101  -1.43701610e-080  -8.33374783e-048   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22430645e-017  -5.75075656e-017   7.99300010e-017]
 [  1.61286945e-022   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  6.04599897e-058   1.69346155e-036   1.53236733e-002   1.25158810e-017]
 [  3.17394363e-105   9.67677792e-084   1.28172714e-049   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429266e-017  -5.75075656e-017  -7.99300010e-017]
 [  2.34172909e-023   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  3.84343394e-060   7.41462611e-038   1.53236733e-002  -1.25158810e-017]
 [ -3.10317155e-109  -6.51628267e-087  -1.97129130e-051   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429066e-017  -5.75075656e-017   7.99300010e-017]
 [  3.39996218e-024   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  2.44326611e-062   3.24640855e-039   1.53236733e-002   1.25158810e-017]
 [  3.03397754e-113   4.38802463e-090   3.03183826e-053   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429037e-017  -5.75075656e-017  -7.99300010e-017]
 [  4.93641339e-025   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  1.55318119e-064   1.42140255e-040   1.53236733e-002  -1.25158810e-017]
 [ -2.96632640e-117  -2.95486877e-093  -4.66295531e-055   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429033e-017  -5.75075656e-017   7.99300010e-017]
 [  7.16719063e-026   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  9.87355328e-067   6.22344715e-042   1.53236733e-002   1.25158810e-017]
 [  2.90018374e-121   1.98979043e-096   7.17160693e-057   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  1.04060616e-026   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  6.27660539e-069   2.72486457e-043   1.53236733e-002  -1.25158810e-017]
 [ -2.83551591e-125  -1.33991263e-099  -1.10299032e-058   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  1.51085862e-027   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  3.99003014e-071   1.19305053e-044   1.53236733e-002   1.25158810e-017]
 [  2.77229004e-129   9.02288913e-103   1.69639478e-060   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  2.19361931e-028   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  2.53645714e-073   5.22363419e-046   1.53236733e-002  -1.25158810e-017]
 [ -2.71047397e-133  -6.07595798e-106  -2.60904850e-062   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  3.18492122e-029   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  1.61242261e-075   2.28710800e-047   1.53236733e-002   1.25158810e-017]
 [  2.65003626e-137   4.09151270e-109   4.01270632e-064   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  4.62419489e-030   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  1.02501502e-077   1.00138386e-048   1.53236733e-002  -1.25158810e-017]
 [ -2.59094618e-141  -2.75519946e-112  -6.17152654e-066   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  6.71387984e-031   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  6.51600760e-080   4.38444378e-050   1.53236733e-002   1.25158810e-017]
 [  2.53317368e-145   1.85533435e-115   9.49178354e-068   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  9.74789851e-032   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  4.14221782e-082   1.91967815e-051   1.53236733e-002  -1.25158810e-017]
 [ -2.47668939e-149  -1.24937073e-118  -1.45983257e-069   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  1.41529976e-032   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  2.63320266e-084   8.40508943e-053   1.53236733e-002   1.25158810e-017]
 [  2.42146457e-153   8.41318556e-122   2.24521676e-071   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  2.05487718e-033   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  1.67392361e-086   3.68007149e-054   1.53236733e-002  -1.25158810e-017]
 [ -2.36747115e-157  -5.66538735e-125  -3.45313454e-073   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  2.98348118e-034   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  1.06411113e-088   1.61127687e-055   1.53236733e-002   1.25158810e-017]
 [  2.31468166e-161   3.81503696e-128   5.31090733e-075   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  4.33172361e-035   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  6.76454108e-091   7.05478996e-057   1.53236733e-002  -1.25158810e-017]
 [ -2.26306927e-165  -2.56902240e-131  -8.16815458e-077   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  6.28924007e-036   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  4.30021024e-093   3.08885843e-058   1.53236733e-002   1.25158810e-017]
 [  2.21260772e-169   1.72996387e-134   1.25625896e-078   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  9.13136299e-037   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  2.73363824e-095   1.35242105e-059   1.53236733e-002  -1.25158810e-017]
 [ -2.16327135e-173  -1.16494702e-137  -1.93212134e-080   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  1.32578482e-037   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  1.73777039e-097   5.92141961e-061   1.53236733e-002   1.25158810e-017]
 [  2.11503508e-177   7.84468147e-141   2.97159502e-082   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  1.92491021e-038   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  1.10469845e-099   2.59262530e-062   1.53236733e-002  -1.25158810e-017]
 [ -2.06787437e-181  -5.28256020e-144  -4.57030145e-084   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  2.79478183e-039   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  7.02255412e-102   1.13515109e-063   1.53236733e-002   1.25158810e-017]
 [  2.02176524e-185   3.55724352e-147   7.02910565e-086   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  4.05775056e-040   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  4.46422881e-104   4.97012816e-065   1.53236733e-002  -1.25158810e-017]
 [ -1.97668425e-189  -2.39542588e-150  -1.08107368e-087   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  5.89145793e-041   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  2.83790463e-106   2.17611330e-066   1.53236733e-002   1.25158810e-017]
 [  1.93260847e-193   1.61306504e-153   1.66268707e-089   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  8.55382214e-042   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  1.80405240e-108   9.52786117e-068   1.53236733e-002  -1.25158810e-017]
 [ -1.88951548e-197  -1.08622807e-156  -2.55720618e-091   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  1.24193152e-042   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  1.14683384e-110   4.17166415e-069   1.53236733e-002   1.25158810e-017]
 [  1.84738337e-201   7.31459291e-160   3.93297304e-093   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  1.80316341e-043   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  7.29040829e-113   1.82651504e-070   1.53236733e-002  -1.25158810e-017]
 [ -1.80619071e-205  -4.92560183e-163  -6.04889706e-095   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  2.61801735e-044   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  4.63450337e-115   7.99718550e-072   1.53236733e-002   1.25158810e-017]
 [  1.76591657e-209   3.31686995e-166   9.30317988e-097   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  3.80110577e-045   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  2.94614796e-117   3.50147546e-073   1.53236733e-002  -1.25158810e-017]
 [ -1.72654045e-213  -2.23355980e-169  -1.43082541e-098   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  5.51883474e-046   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  1.87286255e-119   1.53308065e-074   1.53236733e-002   1.25158810e-017]
 [  1.68804233e-217   1.50406542e-172   2.20060385e-100   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  8.01280961e-047   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  1.19057636e-121   6.71241686e-076   1.53236733e-002  -1.25158810e-017]
 [ -1.65040264e-221  -1.01282839e-175  -3.38452006e-102   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  1.16338178e-047   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  7.56847895e-124   2.93895432e-077   1.53236733e-002   1.25158810e-017]
 [  1.61360223e-225   6.82032402e-179   5.20537853e-104   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  1.68911685e-048   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  4.81127254e-126   1.28678726e-078   1.53236733e-002  -1.25158810e-017]
 [ -1.57762239e-229  -4.59276420e-182  -8.00585170e-106   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  2.45243287e-049   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  3.05851989e-128   5.63404964e-080   1.53236733e-002   1.25158810e-017]
 [  1.54244482e-233   3.09273913e-185   1.23129684e-107   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  3.56069326e-050   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  1.94429725e-130   2.46680367e-081   1.53236733e-002  -1.25158810e-017]
 [ -1.50805164e-237  -2.08263149e-188  -1.89372970e-109   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  5.16977923e-051   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  1.23598731e-132   1.08006154e-082   1.53236733e-002   1.25158810e-017]
 [  1.47442534e-241   1.40243122e-191   2.91254883e-111   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  7.50601507e-052   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  7.85715578e-135   4.72892492e-084   1.53236733e-002  -1.25158810e-017]
 [ -1.44154885e-245  -9.44388550e-195  -4.47948864e-113   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  1.08980016e-052   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  4.99478403e-137   2.07050525e-085   1.53236733e-002   1.25158810e-017]
 [  1.40940543e-249   6.35945435e-198   6.88943590e-115   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  1.58228351e-053   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  3.17517792e-139   9.06546848e-087   1.53236733e-002  -1.25158810e-017]
 [ -1.37797873e-253  -4.28241740e-201  -1.05959253e-116   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  2.29732129e-054   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  2.01845661e-141   3.96921084e-088   1.53236733e-002   1.25158810e-017]
 [  1.34725279e-257   2.88375351e-204   1.62964915e-118   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  3.33548638e-055   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  1.28313032e-143   1.73787320e-089   1.53236733e-002  -1.25158810e-017]
 [ -1.31721196e-261  -1.94190186e-207  -2.50639397e-120   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  4.84280080e-056   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  8.15684331e-146   7.60907740e-091   1.53236733e-002   1.25158810e-017]
 [  1.28784099e-265   1.30766476e-210   3.85482407e-122   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  7.03127430e-057   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  5.18529502e-148   3.33154680e-092   1.53236733e-002  -1.25158810e-017]
 [ -1.25912492e-269  -8.80573395e-214  -5.92870426e-124   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  1.02087243e-057   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  3.29628552e-150   1.45867935e-093   1.53236733e-002   1.25158810e-017]
 [  1.23104915e-273   5.92972702e-217   9.11832382e-126   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  1.48220717e-058   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  2.09544455e-152   6.38665933e-095   1.53236733e-002  -1.25158810e-017]
 [ -1.20359942e-277  -3.99304167e-220  -1.40239462e-127   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  2.15202021e-059   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  1.33207146e-154   2.79632515e-096   1.53236733e-002   1.25158810e-017]
 [  1.17676176e-281   2.68888967e-223   2.15687740e-129   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  3.12452339e-060   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  8.46796153e-157   1.22433872e-097   1.53236733e-002  -1.25158810e-017]
 [ -1.15052251e-285  -1.81068175e-226  -3.31726896e-131   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  4.53650312e-061   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  5.38307248e-159   5.36062591e-099   1.53236733e-002   1.25158810e-017]
 [  1.12486835e-289   1.21930194e-229   5.10194661e-133   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  6.58655995e-062   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  3.42201240e-161   2.34708824e-100   1.53236733e-002  -1.25158810e-017]
 [ -1.09978622e-293  -8.21070420e-233  -7.84677382e-135   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  9.56304245e-063   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  2.17536897e-163   1.02764552e-101   1.53236733e-002   1.25158810e-017]
 [  1.07526337e-297   5.52903765e-236   1.20683073e-136   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  1.38846047e-063   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  1.38287932e-165   4.49942742e-103   1.53236733e-002  -1.25158810e-017]
 [ -1.05128732e-301  -3.72321991e-239  -1.85610091e-138   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  2.01590914e-064   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  8.79094641e-168   1.97002241e-104   1.53236733e-002   1.25158810e-017]
 [  1.02784588e-305   2.50719337e-242   2.85467590e-140   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  2.92690339e-065   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  5.58839357e-170   8.62551598e-106   1.53236733e-002  -1.25158810e-017]
 [ -1.00492714e-309  -1.68832857e-245  -4.39048030e-142   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  4.24957817e-066   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  3.55253475e-172   3.77658271e-107   1.53236733e-002   1.25158810e-017]
 [  9.82519441e-314   1.13691007e-248   6.75254142e-144   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  6.16997291e-067   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  2.25834187e-174   1.65353319e-108   1.53236733e-002  -1.25158810e-017]
 [ -9.60611341e-318  -7.65588239e-252  -1.03853821e-145   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  8.95819872e-068   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  1.43562508e-176   7.23980441e-110   1.53236733e-002   1.25158810e-017]
 [  9.38724727e-322   5.15542406e-255   1.59726769e-147   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  1.30064306e-068   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  9.12625058e-179   3.16986488e-111   1.53236733e-002  -1.25158810e-017]
 [  0.00000000e+000  -3.47163082e-258  -2.45659142e-149   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  1.88840683e-069   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  5.80154600e-181   1.38788879e-112   1.53236733e-002   1.25158810e-017]
 [  0.00000000e+000   2.33777482e-261   3.77822794e-151   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  2.74178247e-070   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  3.68803549e-183   6.07671101e-114   1.53236733e-002  -1.25158810e-017]
 [  0.00000000e+000  -1.57424317e-264  -5.81089971e-153   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  3.98080066e-071   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  2.34447952e-185   2.66061785e-115   1.53236733e-002   1.25158810e-017]
 [  0.00000000e+000   1.06008566e-267   8.93714089e-155   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  5.77973419e-072   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  1.49038268e-187   1.16492084e-116   1.53236733e-002  -1.25158810e-017]
 [  0.00000000e+000  -7.13855155e-271  -1.37452875e-156   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  8.39161016e-073   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  9.47434391e-190   5.10047159e-118   1.53236733e-002   1.25158810e-017]
 [  0.00000000e+000   4.80705666e-274   2.11401981e-158   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  1.21837993e-073   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  6.02282850e-192   2.23318267e-119   1.53236733e-002  -1.25158810e-017]
 [  0.00000000e+000  -3.23704236e-277  -3.25135416e-160   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  1.76896879e-074   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  3.82870449e-194   9.77773282e-121   1.53236733e-002   1.25158810e-017]
 [  0.00000000e+000   2.17980439e-280   5.00056990e-162   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  2.56837009e-075   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  2.43390263e-196   4.28106757e-122   1.53236733e-002  -1.25158810e-017]
 [  0.00000000e+000  -1.46786686e-283  -7.69085684e-164   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  3.72902280e-076   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  1.54722883e-198   1.87441607e-123   1.53236733e-002   1.25158810e-017]
 [  0.00000000e+000   9.88452503e-287   1.18285076e-165   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  5.41417730e-077   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  9.83571422e-201   8.20691462e-125   1.53236733e-002  -1.25158810e-017]
 [  0.00000000e+000  -6.65617829e-290  -1.81921981e-167   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  7.86085723e-078   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  6.25255115e-203   3.59330293e-126   1.53236733e-002   1.25158810e-017]
 [  0.00000000e+000   4.48222948e-293   2.79795291e-169   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  1.14131978e-078   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  3.97473889e-205   1.57328626e-127   1.53236733e-002  -1.25158810e-017]
 [  0.00000000e+000  -3.01830573e-296  -4.30324057e-171   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  1.65708497e-079   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  2.52673651e-207   6.88845246e-129   1.53236733e-002   1.25158810e-017]
 [  0.00000000e+000   2.03250850e-299   6.61836707e-173   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  2.40592571e-080   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  1.60624322e-209   3.01602948e-130   1.53236733e-002  -1.25158810e-017]
 [  0.00000000e+000  -1.36867871e-302  -1.01790225e-174   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  3.49316940e-081   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  1.02108679e-211   1.32053373e-131   1.53236733e-002   1.25158810e-017]
 [  0.00000000e+000   9.21659814e-306   1.56552966e-176   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  5.07174117e-082   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  6.49103587e-214   5.78180466e-133   1.53236733e-002  -1.25158810e-017]
 [  0.00000000e+000  -6.20640044e-309  -2.40777845e-178   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  7.36367338e-083   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  4.12634333e-216   2.53149650e-134   1.53236733e-002   1.25158810e-017]
 [  0.00000000e+000   4.17935184e-312   3.70315379e-180   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]
 [  1.06913353e-083   3.49984625e-001   5.16395464e-017   1.76827625e-017]
 [  2.62311126e-218   1.10838655e-135   1.53236733e-002  -1.25158810e-017]
 [  0.00000000e+000  -2.81434980e-315  -5.69543596e-182   2.35677492e-004]]
current a:
[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]
 [  1.55227759e-084   3.49984625e-001   5.16395464e-017  -1.76827625e-017]
 [  1.66750853e-220   4.85294271e-137   1.53236733e-002   1.25158810e-017]
 [  0.00000000e+000   1.89516665e-318   8.75955810e-184   2.35677492e-004]]

Convergence and complexity of the QR-algorithm

The convergence of the QR-algorithm is from largest eigenvalues to the smallest.

At least 2-3 iterations is needed for an eigenvalue.

Each step is one QR-factorization and one matrix-by-matrix product. $\mathcal{O}(n^3)$ complexity.

It means, $\mathcal{O}(n^4)$ complexity?

Fortunately, not.

We can also speedup the QR-algorithm by using shifts.

A few words about the SVD

And the last, but not the least, is the singular value decomposition of matrices.

$$A = U \Sigma V^*.$$

We can compute via eigendecomposition of

$$A^* A = U^* \Sigma^2 U,$$

but it is a bad idea (c.f. Gram matrices)

We will discuss methods for computing SVD later (there will a special SVD lecture).

Summary of todays lecture

  • LU decomposition and Gaussian elimination
  • QR decomposition and Gram-Schmidt algorithm, reduction to a simpler form by Householder transformations
  • Schur decomposition and QR-algorithm

Next lecture

  • No Test
  • SVD: how to compute it, what are the applications, computing best rank-$r$ approximation
Questions?
In [42]:
from IPython.core.display import HTML
def css_styling():
    styles = open("./styles/custom.css", "r").read()
    return HTML(styles)
css_styling()
Out[42]:
In [ ]: