Square systems

Singularity, Degenerate Systems

If the rows are linear combinations of one another, or the columns linear combinations or one another then even if there seems to be the same number of unknowns as equations, the system is in fact equivalent to a singular set of equations.

Close singularity is just as bad because the roundoff could take us there, we will have to measure this.

Single precision is usually enough for 30-50 equations.

Double precision makes it possible to solve for a few hundred.

This does NOT hold for pathological systems we will look into these and how to recognize them when we talk about conditionning.

A full square system

:

Suppose the matrix of coefficients is a square matrix, (of full rank also, as we saw, non-degenerate) then the equation looks like this:

$$Ax=b$$

Formally you would say : invert it, in fact a computational mistake, except if you have multiple right hand sides.

Doing what one would first think of doing if you had to solve a system by hand.

The easiest system: triangular.

Bringing to triangular shape.

Pivotting. This is NOT to solve ill-conditionning, but unlucky diagonal elements that are very small make the algorithm unstable.

LU decomposition, with a permutation matrix (prepivotting) is what is best to do.

$$A=LU, L\; lower\;triang.\;U\; upper\; triang.$$

Flops: Of an order $\frac{n^3}{2}$ (the same as the cost of the multiplication of two matrices, not so bad).

Condition number for the inversion problem

Try a small perturbation of $b$ to $b^*=b+\delta b$, call the solution of the perturbed system: $x^*=x+\delta x$,

$$\begin{eqnarray*} A(x+\delta x)&=&b+\delta b\\ x+\delta x&=&A^{-1}(b+\delta b)\... ...vert \frac{\vert\vert\delta b\vert\vert}{\vert\vert b\vert\vert} \end{eqnarray*}$$

The condition number for the data $A$ and the linear solving problem is thus:

$$\kappa=\vert\vert A \vert\vert \vert\vert A^{-1} \vert\vert$$

(depends on the norm, not on the algorithm). For the two norm we have a better explanation of how to compute this, we need the following to do it:

Singular Value Decomposition

This is the most important decomposition in statistics.

We will study it several times, here is my first introduction: a little made-up example with matlab to start with:

>> u=[3 1 -1 2]'
u =
     3
     1
    -1
     2
>> v=(1:4)'
v = 
     1
     2
     3
     4
>> X=u*v'
X =
     3     6     9    12
     1     2     3     4
    -1    -2    -3    -4
     2     4     6     8
>>svd(X)
ans =
   21.2132
    0.0000
         0
         0

>> E=10^(-3)*randn(4);
>> XE=X+E 
XE =
 
    3.0012    5.9993    9.0003   12.0012
    1.0006    2.0017    3.0009    3.9994
   -0.9999   -1.9999   -3.0014   -3.9994
    2.0004    4.0018    5.9993    7.9996

>> svd(XE) 
ans =
   21.2143
    0.0028
    0.0019
    0.0005
>> cond(X)
Condition is infinite
ans =
   Inf 
>> cond(XE)
ans =
 
   4.2567e+04
%An example you can't see with your bare eyes:
>> X2=u'*v 
X2 =
    0.2468    0.2499    0.0166    0.2487
    0.3225    0.3266    0.0218    0.3250
    0.1495    0.1514    0.0101    0.1506
    0.4367    0.4423    0.0295    0.4401
 
>> svd(X2) 
ans =
 
    1.0730
    0.0000
         0
         0
>> [U,S,V]=svd(X2);
>> U*S*V' 
ans =
    0.2468    0.2499    0.0166    0.2487
    0.3225    0.3266    0.0218    0.3250
    0.1495    0.1514    0.0101    0.1506
    0.4367    0.4423    0.0295    0.4401
>> 10000*(X2-U*S*V') 
ans =
    1.0e-11 *
   -0.0278    0.0833   -0.0035    0.0555
   -0.0555    0.0555         0    0.0555
         0    0.0555    0.0017    0.0278
   -0.0555    0.1665         0    0.1110

>> A=rand(4)
A = 
    0.5045    0.4940    0.0737    0.9138
    0.5163    0.2661    0.5007    0.5297
    0.3190    0.0907    0.3841    0.4644
    0.9866    0.9478    0.2771    0.9410
>> flops(0)
>> [L,U,P]=lu(A) 
L =
     1.0000         0         0         0
    0.5233    1.0000         0         0
    0.5114   -0.0406    1.0000         0
    0.3234    0.9388    0.7363    1.0000
U =
     0.9866    0.9478    0.2771    0.9410
         0   -0.2298    0.3557    0.0373
         0         0   -0.0535    0.4342
         0         0         0   -0.1945
>> flops 
ans =
    34 
>> P*A
ans =
    0.9866    0.9478    0.2771    0.9410
    0.5163    0.2661    0.5007    0.5297
    0.5045    0.4940    0.0737    0.9138
    0.3190    0.0907    0.3841    0.4644
>> L*U
ans =
     0.9866    0.9478    0.2771    0.9410
    0.5163    0.2661    0.5007    0.5297
    0.5045    0.4940    0.0737    0.9138
    0.3190    0.0907    0.3841    0.4644
>> P*A-L*U
ans =    1.0e-15 *
         0         0         0         0
   -0.1110         0         0         0
         0         0         0         0
   -0.0555   -0.0139         0         0

Here is what we need to remember:

$$X=USV', V'V=I, U'U=I, S\; diagonal\; s_i$$

Actually the singular values are the square roots of the eigenvalues of $X'X$.

The condition number in $\ell^2$ norm has the simple interpretation as a function of the singular values, ration of the largest to the smallest:

$$\kappa=\frac{s_1}{s_n}$$

Read 2.6 of NR (pages 61-71).

Tridiagonal

Of an order of $N$ steps are sufficient, see NR 51.

Iterative Improvement

This enables, a the price of storing A, to improve the solution by using the residual. Do the following in double precision:

$$\begin{eqnarray*} r=b-Ax^*,solve \;\; Ly&=&Pr\;\; for\;y\; \\ and\; \; Uz&=&y\; for\; z \\ x_{new} &= &x^*+z\\ Ax_{new}&=&A(x^*+z)=(b-r)+Az=b \end{eqnarray*}$$

Cholesky Decomposition

For a symmetric positive definite matrix A the LU decomposition can be rewritten:

$$A=GG'=(LD^{1/2})(D^{1/2}L')$$

with D a diagonal matrix with only positive entries.