C<sup>4</sup>E: Computation For Economists

Linear Systems

Love Triangles

$A = \l(\matrix{ 1 & 0 & 0 \cr 2 & 3 &0\cr 4 & 5 & 6\cr}\r).\nonumber$

$\l(\matrix{ 1 & 0 & 0 \cr 2 & 3 &0\cr 4 & 5 & 6\cr}\r)\l(\matrix{x\cr y\cr z}\r) = \l(\matrix{-1 \cr 4 \cr 0}\r).\nonumber$

$x$

$x = -1/1 = -1$

$x$

$y$

$2(-1) + 3y = 4$

$y = 6/3 = 2$

$z$

$4(-1)+5(2) + 6z = 0$

$z = (0+4-10)/6 = -1.$

$N$

$a_{ii}$

$A$

$A \hbox{ triangular } \quad \rightarrow\quad |A| = \prod_{n=1}^N a_{nn}.\nonumber$

computing

$a_{ii}=0$

$a_{ii}$

overflow

feq(A[n][n],0.0)

Algorithm 3. Solving a Lower Triangular Systems.

Suppose

$A$ is a

$N\times N\$ lower triangular matrix and we wish to solve

$Ax = b\nonumber$ for

$x$ .

This requires a sequence of forward substitutions

$\matrix{& x_1 &= b_1/a_{11} \cr \downarrow&x_2 &= {b_2-a_{21}x_1 \over a_{22}} \cr \downarrow & \vdots &\vdots \cr \downarrow& x_N &= {b_N-S_N \over a_{NN} }\cr}\nonumber$

At row

$n$ , the term

$S_n = \sum_{i=1}^{n-1} a_{ni}x_{i}\nonumber$ is the sum of values of earlier rows substituted into the current

$n^{th}$ equation. When

$n=1$ the sum is simply

$S_1=0$ .

The elements of

$x$ are computed by looping over rows of

$A$ , and at each row the sum

$S_n$ loops over columns of the current row up to

$n$ .

nested loop

$D=1$

$n$

$a_{nn}$

numerically

$D *= a_{nn}$

$x_{n}$

$0$

$D=0$

$N$

$N(N-1)$

$N+N(N-1) = N^{2}$

FLOPS

$A$

upper triangular

$x$

back substitutions

$^{2}$

$A$

$2N^{2}$

$x_N = b_{N} / a_{NN}$

$x_N$

$N-1$

$n=1$

Algorithm 4. Solving Upper Triangular Systems.

Suppose

$A$ is a

$N\times N\$ upper triangular matrix and we wish to solve

$Ax = b\nonumber$ for

$x$ .

This requires a sequence of back substitutions

$\matrix{a_{11}&a_{12}&\ldots&a_{1n}\cr 0&a_{22}&a_{23}&\ldots\cr 0&0&\ddots&\ldots\cr 0&\ldots&0&a_{nn}\cr}x = b \quad\Rightarrow\quad \matrix{& x_1 &=& {b_1-\sum_{i=2}^{n} a_{1i}x_{i}\over a_{11} }\cr \uparrow &\vdots \cr \uparrow& x_{n-1} &=&{b_{n-1}-a_{n-1 n}x_n \over a_{n-1 n-1}} \cr \uparrow& x_n &=& b_n/a_{nn}. \cr }\nonumber$

Beethoven in his grave … Decomposing

So we have simple algorithms to solve any triangular system. Now Suppose $A = LU$ , where $L$ is lower $\triangle$ and $U$ is upper $\triangle$ . Then $Ly = b \quad\Rightarrow\quad y = L^{-1}b \qquad Ux = y \quad\Rightarrow\quad x = U^{-1}y = U^{-1}L^{-1}b\nonumber$ Now, to consider the general problem consider first the trivial case when $N=1$ (one equation in one unknown). Then we have scalars $Ax = b$ and therefore $x = b/A$ . It should be clear that any non-zero scalar A can be written as the product of two non-zero numbers, call them $L$ and $U$ . Indeed, in the scalar case we can simply normalize $L=1$ and define $U=A$ . When $N=1$ there is a decomposition of $A$ into two triangular "matrices". If, on the other hand, $A=0$ , then $U=0$ and it is not possible to factor the matrix and then solve, since this would involve dividing by 0.

Once we let $N$ be greater than 1 it is not so obvious that $A$ can be factored into triangular matrices or how to go about doing it. However, the idea that elements of $L$ can still be normalized might be intuitively clear. There are $N(N+1)/2$ non-zero elements in a triangular matrix. Two such matrices have a total of $N(N+1) = N^2+N$ elements, but the original $A$ matrix has only $NN$ elements. So we have $N$ more degrees of freedom (free values) than we need for the factorization. One way to reduce this freedom is to normalize $N$ elements of the two triangles. Following the scalar example, we can normalize the diagonal of $L$ to consist of 1s. Then together there are exactly $N^2$ free elements in L and U to just match the elements in A.

Fortunately, for the case of $N>1$ it has been known for decades that if $A$ is non-singular then we can indeed find matrices $L$ and $U$ for it. We will not discuss how the algorithm for computing the $L$ and $U$ works (it's tedious) but let's be clear what it does.

Algorithm 5. L-U Decomposition & solving $Ax=b$ & Inversion

NR 2.3; Judd II.3

$|A|\ne 0$ then:

$\exists L, U: A=LU$ and $L$ is lower $\triangle$ , $U$ is upper $\triangle$ , and $L$ has 1s on its diagonal.
Since the diagonal of $L$ is normalized $|A| = |LU| = |L| |U| = 1\times |U| = \prod_{i=1,\dots N} u_{ii}.$
The matrices $L$ and $U$ for $A$ can be found numerically using an algorithm that fails only if $|A|=0$ . This will happen in the process at some row $i$ when $u_{ii}$ is too close to 0. At that point the L-U decomposition will stop because it has discovered that the determinant of $A$ is 0. Thus, instead of computing the determinant in order to solve a system of equations (as taught in math econ courses), computer algorithms compute the determinant as a by-product of solving the system.
The solution is numerically more accurate if permutations of the rows of $A$ are allowed. That is, a permutation matrix $P$ is the identity matrix $I$ with some of it is rows swapped. The L-U algorithm will then produce 3 matrices such that $A = PLU$ . However, note Ox's implements this differently

$A$

B

B[<1;0>][]

swaps

The built-in Ox function for carrying out the LU Decomposition is:

declu(A, aL, aU, aP)

        Argument        Explanation
        ------------------------------------------------------------------------------------
        A               the (square) matrix to decompose
        aL              an address (pointer) where to put L
        aU              an address (pointer) where to put U
        aP              an address (pointer) where to put a permutation matrix P,

        On output
        P contains a 2xN matrix containing integers.  The first row is the order of rows of LU to get back to the original rows of A.
            That is, A = (L*U)[ P[0][] ][] 

  Returns
    1:  if no error occurred, and the L, U, P matrices are reliable.
    2:  if the decomposition could be unreliable (because A is near to singular)
    0:  if the LU decomposition failed and the matrix A is (numerically) singular.

declu()

declu(A,&L,&U,&P)

$Ax = b$

solvelu()

    decl L, U, P, ok, A, b;

    // compute or load A and b

    ok = declu(A, &L, &U, &P);
    if (ok>0) {
        x = solvelu(L,U,P,b);
        println(" x = ",x);
        }
    else println("A is singular, cannot solve for x");

One key implication of solving linear systems using decompositions is that it does not require computing $A^{-1}$ , although if $A^{-1}$ is already available it can be used to do the job. Since $L-U$ depends only on $A$ , the decomposition can be used to solve for different $b$ 's. In fact, suppose you 5 different $b$ vectors to solve $Ax=b$ . Then you can put the column vectors into a matrix $B$ . The back/forward substitution process can be used to solve a matrix equation $AX=B$ . Each column of $X$ will correspond to the vector solution for the corresponding column of $B$ .

$A$

$N\times N$

$B = I$

$N$

$A^{-1}$

$AX = I$

$A^{-1}$

$n-1$

$N$

$A^{-1}$

$b$

$N$

$A^{-1}$

$A^{-1}b$

$A^{-1}$

$b$

declu()

solvelu()

invert(A)

    Returns
        inverse of A,
        0 if decomposition failed.

invert(A,aLogDet,aSign)

        Argument     Explanation
        ------------------------------------------------------------------------------------
        aLogDet      address to put log( |det(A)| ), the logarithm of absolute value of |A|
        aSign        address to put a sign code for |A|:
                            0: singular
                        -1,-2: negative determinant
                        +1,+2: positive determinant

                        either -2 or +2 indicate the result is unreliable.

Algorithm 6. A General Algorithm to Solve $Ax=b$ & Inversion of a Matrix

Use the L-U decomposition to compute $P$ , $L$ and $U$ for $A$ . If this fails then $A$ is not invertible. or handle this situation some other way.
Forward solve $Ly=b$ for $y$ .
Backward solve $Uz = y$ for $z$ .
Rearrange rows: $x = Pz$ .

Non-square and Singular Systems: Singular Value Decomposition

$B$

$N\times M$

$N\ge M$

$N\ne M$

$B$

$B = UWV'.\tag{SD}\label{SVD}$

$U$

$N\times M$

$B$

$W$

$M\times M$

$V$

$M\times M$

$W$

$M$

$w_i.$

singular values

$B.$

$U$

$V$

$U'U = V'V = I_M.$

$U$

$V$

$B$

$N=M$

$B^{-1} = V\hbox{diag}\pmatrix{1/w_1 & 1/w_2 & \cdots & 1/w_N}U'.$

$UU'=I$

$B$

$w's$

$N>M$

$Bx = b.$

$x$

$x = V\hbox{diag}\pmatrix{1/w_1 & 1/w_2 & \cdots & 1/w_N}U'b.$

$x = B^{-1}b.$

$x$

$Bx$

$b.$

$b \equiv y,$

$\hat\beta^{ols} \equiv x$

$X$

$(X'X)^{-1}X'.$

Exercises

Write a nested-loop flowchart for solving a triangular system. The outer loop is over rows; the inner loop is over columns of the current row.
Write a function that solves a triangular system. Test your function by sending triangular systems that you can verify easily. (In this version just assume the matrix is non-singular.) You function should satisfy this description:

SolveL(A,b,ax)

        Argument        Explanation
        ------------------------------------------------------------------------------------
        A               N × N lower triangular matrix
        b               N × 1 constant vector
        ax              address to place resulting solution x*

23-SolveL.ox
 1:    #include 
 2:    
 3:    SolveL(A,b,ax) {
 4:    	 decl n, j, S;
 5:    	 ax[0]  = constant(.NaN,b); 	  //initialize x* 		
 6:    	 for(n=0;n < rows(A);++n) {
 7:             // put code here to solve ax[0][n]
 8:    		}
 9:    	}
10:    
11:    main() {
12:    	// Put calls to your function to test it
13:    
14:    	}

Once your SolveL() function works, modify it so that it implements an additional feature:

⋮
        Returns
        the determinant of A, |A|.  If A has a numerically zero determinant it returns exactly 0 as the determinant and does not complete the solution to x*.

Once your SolveL() code works try eliminating the inner loop by using vector multiplication and subranges of A and ax[0].
Make a copy of your SolveL() function, rename it SolveU(), and change the code to solve an upper triangular system.
Write your own MyInvert(A) function to replace the standard invert(A) using declu() and solvelu(). Your function returns the inverse if successful or the integer 0 if the inverse fails.
Modify MyInvert(A) to not use solvelu(). Instead, use your SolveL() and SolveU() coded earlier.

Linear Systems

Algorithm 3. Solving a Lower Triangular Systems.

Algorithm 4. Solving Upper Triangular Systems.

Algorithm 5. L-U Decomposition & solving Ax=bAx=b & Inversion

Algorithm 6. A General Algorithm to Solve Ax=bAx=b & Inversion of a Matrix

Exercises

Algorithm 5. L-U Decomposition & solving $Ax=b$ & Inversion

Algorithm 6. A General Algorithm to Solve $Ax=b$ & Inversion of a Matrix