<jupyter_start><jupyter_text>2 Linear Equations and Computer Basics<jupyter_code>#%pylab inline
%pylab notebook
# pylab Populating the interactive namespace from numpy and matplotlib
# numpy for numerical computation
# matplotlib for ploting
import numpy as np
from numpy import append, array, diagonal, tril, triu
from numpy.linalg import inv
from numpy.linalg import norm
from numpy.linalg import solve
from scipy.linalg import lu
#from scipy.linalg import solve
from pprint import pprint
from numpy import array, zeros, diag, diagflat, dot

from sympy import *
import sympy as sym
init_printing()<jupyter_output><empty_output><jupyter_text>2.1 Linear Equations and the L-U FactorizationThe linear equation is the most elementary problem that arises in computationaleconomic analysis. In a linear equation, an $n \times n$ matrix A and an n-vector b aregiven, and one must compute the n-vector x that satisfies$$A x = b$$ Most linear equations encountered in practice, however, do not have a triangular A matrix. In such cases, the linear equation is often best solved using the L-U factorization algorithm. The L-U algorithm is designed to decompose the A matrix into the product of lower and upper triangular matrices, allowing the linear equation to be solved using a combination of backward and forward substitution. The L-U factorization algorithm is faster than other linear equation solution methods that are typically presented in elementary linear algebra courses.<jupyter_code># example problem

A = np.array([[-3.0, 2.0, 3.0], [-3.0, 2.0, 1.0], [3.0, 0, 0]])

b = [10, 8, -3]

x0 = [1.0, 1.0, 1.0]

n = 25


print( "A:")
pprint(A)


#rint( "diag(A):")
#print((diag(A)))


#rint("pprint(diagflat(diag(A)))")
#print(diagflat(diag(A)))

print( "b:")
pprint(b)

print( "x0:")
pprint(x0)

print( "n:")
pprint(n)<jupyter_output>A:
 [[-3.  2.  3.]
  [-3.  2.  1.]
 [ 3.  0.  0.]]
b:
[10, 8, -3]
x0:
[1.0, 1.0, 1.0]
n:
25<jupyter_text>Using Sympy<jupyter_code>M = Matrix(A)
v = Matrix(b)
M
v<jupyter_output><empty_output><jupyter_text>http://docs.sympy.org/latest/modules/matrices/matrices.htmlReturns (L, U, perm) where L is a lower triangular matrix with unit diagonal, U is an upper triangular matrix, and perm is a list of row swap index pairs. If A is the original matrix, then A = (L*U).permuteBkwd(perm), and the row permutation matrix P such that P*A = L*U can be computed by P=eye(A.row).permuteFwd(perm).See documentation for LUCombined for details about the keyword argument rankcheck, iszerofunc, and simpfunc.<jupyter_code>M.LUdecomposition()
L,U,perm = M.LUdecomposition()


L
U
perm
L*U
(L*U).permuteBkwd(perm)
v = v.permuteBkwd(perm)
v<jupyter_output><empty_output><jupyter_text>the elements of x can be computed recursively using forward-substitution<jupyter_code>L.LUsolve(v)<jupyter_output><empty_output><jupyter_text>If A is an upper triangular matrix, then the elements of x can be computed recursively using backward-substitution.<jupyter_code>U.LUsolve(L.LUsolve(v))
y= b
b = [10, -3, 8]
y[0] = b[0]/L[0,0]
for i in range(len(b)-1):
    y[i+1] = (b[i+1] - L[i+1,:i+1].dot(y[:i+1]))/L[i+1,i+1]  
    
y    
L.dot([10, 7.0 , -2.0])<jupyter_output><empty_output><jupyter_text>Using Numpy<jupyter_code>"""
Solve a linear equation by LU-decomposition
Comes from LU decomposition of a matrix A s.t. A = LU
Then
LUx = b => Ux = y => Ly = b
"""
p, l, u = lu(A)
p

l@u
p@l@u<jupyter_output><empty_output><jupyter_text>Lower triangle matrix<jupyter_code>l
p@l
b=p@l.dot([10, 7, -2])
b<jupyter_output><empty_output><jupyter_text>Forward solving<jupyter_code>y = np.linalg.solve(p@l, b)
y<jupyter_output><empty_output><jupyter_text>Upper triangle matrix Backward solving<jupyter_code>x = np.linalg.solve(u, y)
x   
## Other example

A = np.array([[54, 14, -11, 2],
               [14, 50, -4, 29],
               [-11, -4, 55, 22],
               [2, 29, 22, 95]])
print("Left hand side of linear system, matrix A:\n", A)
b = np.array([1, 1, 1, 1])
x0 = [1, 1, 1,1]

n = 20

p, l, u = lu(A)

print("perm matrix:\n",p)

print("Lower matrix:\n",l)

print("Upper matrix\n",u)

print("Right hand side of linear equation system:\n",b)

y = np.linalg.solve(l, b)

print("Intermedia output matrix y:\n",y)    


x = np.linalg.solve(u, y)
print("Solution x:\n",x)<jupyter_output>Left hand side of linear system, matrix A:
 [[ 54  14 -11   2]
 [ 14  50  -4  29]
 [-11  -4  55  22]
 [  2  29  22  95]]
perm matrix:
 [[ 1.  0.  0.  0.]
 [ 0.  1.  0.  0.]
 [ 0.  0.  1.  0.]
 [ 0.  0.  0.  1.]]
Lower matrix:
 [[ 1.          0.          0.          0.        ]
 [ 0.25925926  1.          0.          0.        ]
 [-0.2037037  -0.02476038  1.          0.        ]
 [ 0.03703704  0.61421725  0.43831321  1.        ]]
Upper matrix
 [[ 54.          14.         -11.           2.        ]
 [  0.          46.37037037  -1.14814815  28.48148148]
 [  0.           0.          52.73083067  23.11261981]
 [  0.           0.           0.          67.30154198]]
Right hand side of linear equation system:
 [1 1 1 1]
Intermedia output matrix y:
 [ 1.          0.74074074  1.22204473 -0.02765113]
Solution x:
 [ 0.01893441  0.01680508  0.02335523 -0.00041085]<jupyter_text>2.2 Gaussian EliminationThe L-U factors of a matrix A are computed using Gaussian elimination. Gaussianelimination is based on two elementary row operations: subtracting a constant multipleof one row of a linear equation from another row, and interchanging two rows ofa linear equation. Either operation may be performed on a linear equation withoutaltering its solution.The Gaussian elimination algorithm begins with matrices L and U initialized asL = I and U = A, where I is the identity matrix. The algorithm then uses elementaryrow operations to transform U into an upper triangular matrix, while preserving thepermuted lower diagonality of L and the factorization A = LU:  2.3 Rounding ErrorIn practice, Gaussian elimination performed on a computer can sometimes renderinaccurate solutions due to rounding errors. The eects of rounding errors, however,can often be controlled by pivoting.  2.4 Ill ConditioningPivoting cannot cure all the problems caused by rounding error. Some linear equationsare inherently diÆcult to solve accurately on a computer, despite pivoting. Thisoccurs when the A matrix is structured in such a way that a small perturbation Æbin the data vector b induces a large change Æx in the solution vector x. In such casesthe linear equation or, more generally, the A matrix are said to be **ill-conditioned**.  2.5 Special Linear EquationsGaussian elimination can be accelerated for matrices possessing certain special structures.Two such classes arising frequently in computational economic analysis are**symmetric positive deffinite matrices and sparse matrices**.  Cholesky factorization algorithmA special formof Gaussian elimination, the Cholesky factorization algorithm, may be applied tosuch linear equations. Cholesky factorization requires only half as many operationsas general Gaussian elimination and has the added advantage that it is less vulnerableto rounding error and does not require pivoting.The essential idea underlying Cholesky factorization is that any symmetric positivedeffinite matrix A can be uniquely expressed as the product$$\mathbf {A=U^T U}$$ $${\displaystyle {\begin{aligned}\left({\begin{array}{*{3}{r}}4&12&-16\\12&37&-43\\-16&-43&98\\\end{array}}\right)=\left({\begin{array}{*{3}{r}}2&0&0\\6&1&0\\-8&5&3\\\end{array}}\right)\left({\begin{array}{*{3}{r}}2&6&-8\\0&1&5\\0&0&3\\\end{array}}\right)\end{aligned}}}$$ $$\mathbf {A s=U^T U x = U^T (U x)     = b}$$may be solved eÆciently by using forward substitution to solve$$\mathbf { U^T y = b}$$and then using backward substitution to solve$$\mathbf { U x = y}$$ https://en.wikipedia.org/wiki/Cholesky_decompositionIn linear algebra, the Cholesky decomposition or Cholesky factorization is a decomposition of a Hermitian, positive-definite matrix into the product of a lower triangular matrix and its conjugate transpose, which is useful e.g. for efficient numerical solutions and Monte Carlo simulations. It was discovered by André-Louis Cholesky for real matrices. When it is applicable, the Cholesky decomposition is roughly twice as efficient as the LU decomposition for solving systems of linear equations.  2.6 Iterative MethodsAlgorithms based on Gaussian elimination are called **exact** or, more properly, **direct**methods because they would generate exact solutions for the linear equation $Ax = b$after a finite number of operations, if not for rounding error. Such methods are idealfor moderately-sized linear equations, but may be impractical for large ones. Othermethods, called **iterative methods** can often be used to solve large linear equationsmore efficiently if the A matrix is sparse, that is, if $A$ is composed mostly of zeroentries. Iterative methods are designed to generate a sequence of increasingly accurateapproximations to the solution of a linear equation, but generally do not yield an exactsolution after a prescribed number of steps, even in theory. The most widely-used iterative methods for solving a linear equation $Ax = b$ aredeveloped by choosing an easily invertible matrix $Q$ and writing the linear equationin the equivalent form$$Q x = b + (Q - A)x$$or$$x = Q^{-1}b + (I - Q^{-1}A)x.$$This form of the linear equation suggests the iteration rule$$x^{(k+1)} \leftarrow   Q^{-1}b + (I - Q^{-1}A)x^{(k)}$$ Ideally, the so-called splitting matrix Q will satisfy two criteria. First, $Q^{-1}b$ and$Q^{-1}A$ should be relatively easy to compute. This is true if Q is either diagonal ortriangular. Second, the iterates should converge quickly to the true solution of thelinear equation. If$$||I - Q^{-1}A|| < 1$$in any matrix norm, then the iteration rule is a contraction mapping and is guaranteedto converge to the solution of the linear equation from any initial value. The two most popular iterative methods are the **Gauss-Jacobi** and **Gauss-Seidel**methods.  - The Gauss-Jacobi method sets $Q$ equal to the diagonal matrix formed from the diagonal entries of A.  - The Gauss-Seidel method sets $Q$ equal to the upper triangular matrix formed from the upper triagonal elements of A.  Jacobi-MethodHere, the user specifies the data A and b and an initial guess $x$ for the solution ofthe linear equation, typically the zero vector or $b$. Iteration continues until the normof the change dx in the iterate falls below the specified convergence tolerance tol oruntil a specified maximum number of allowable iterations maxit are performed.<jupyter_code>#p17

maxit = 1000
tol = 1/10000 
x0 = np.array([1,1,1,1])
x0= b

d = np.diag(np.diag(A))
x = x0
for it in np.arange(maxit):
    dx =  inv(d).dot(b - A.dot(x))
    x = x + dx
    if np.linalg.norm(dx)<tol:
        break

print("x:\n", x)        
def gjacobi(A, b, maxit=1000,tol = 1/1000, x0=None):
    """
    Solve a linear equation by the gauss jacobi iteration outlined in the book.
    Follows the eq:
        x = inv(D)(b - Rx)
    Where D is the diagonal matrix of A and R is the remainder s.t D + R = A
    """
    # If we have not provided an initial array for x make a new one
    if x0==None:
        x = array([1 for _ in range(A.shape[1])])
    else:
        x = x0
    d = np.diag(np.diag(A))
    
    for it in np.arange(maxit):
        dx =  inv(d).dot(b - A.dot(x))
        x = x + dx
        if np.linalg.norm(dx)<tol:
            break

    print("x:\n", x)  
gjacobi(A, b, maxit,tol, x0=None)<jupyter_output>x:
 [ 0.01892514  0.01678369  0.02335003 -0.00042463]<jupyter_text>Jacobi-Method-1The Jacobi method is a matrix iterative method used to solve the equation $Ax=b$ for a known square matrix A of size n×n and known vector b or length n.$$\begin{eqnarray*}Ax = b\end{eqnarray*}$$A is split into the sum of two separate matrices, $D$ and $R$, such that $A=D+R$. $D_{ii}=A_{ii}$, but $D_{ij}=0$, for $i≠j$. $R$ is essentially the opposite. $R_{ii}=0$, but $R-{ij}=A_{ij}$ for $i≠j$. The solution to the equation, i.e. the value of $x$, is given by the following iterative equation:$$\begin{eqnarray*}x^{(k+1)} = D^{-1}(b-Rx^{(k)}).\end{eqnarray*}$$https://www.quantstart.com/articles/Jacobi-Method-in-Python-and-NumPy<jupyter_code>def gauss_jacobi(A, b, maxit, x0=None):
    """
    Solve a linear equation by the gauss jacobi iteration outlined in the book.
    Follows the eq:
        x = inv(D)(b - Rx)
    Where D is the diagonal matrix of A and R is the remainder s.t D + R = A
    """
    D = np.diag(np.diag(A))
    # Calculate the remainder matrix
    R = A - D
    # If we have not provided an initial array for x make a new one
    if x0==None:
        x = array([1 for _ in range(A.shape[1])])
    else:
        x = x0
    for _ in range(maxit):
        x = inv(d).dot(b - R.dot(x))
    return x

gauss_jacobi(A, b, maxit, x0)<jupyter_output>/opt/conda/lib/python3.5/site-packages/ipykernel/__main__.py:12: FutureWarning: comparison to `None` will result in an elementwise object comparison in the future.<jupyter_text>Check result with direct method<jupyter_code>np.linalg.solve(A,b)<jupyter_output><empty_output><jupyter_text>gauss seidel methodhttps://docs.scipy.org/doc/numpy-dev/user/numpy-for-matlab-users.html<jupyter_code>A = np.array([[3, 1],
               [2, 5]])
b = np.array([7 , 9])
Q = np.tril(A)

x=b
lamb = 1.5

for it in np.arange(maxit):
    dx = np.linalg.inv(Q).dot(b -A.dot(x))
    x = x + lamb*dx
    if np.linalg.norm(dx)<tol:
        break
print("x:\n", x)     

def gseidel(A, b, maxit=1000, tol=1/1000,lamb =1, x0=None):
    """
    Solve a linear equation by the gauss seidel iteration outlined in the book
    Follows the eq:
        
    """
    Q = np.tril(A)
    # If we have not provided an initial array for x make a new one
    if x0 == None:
        x = array([1 for _ in range(A.shape[1])])
    else:
        x = x0
    for it in np.arange(maxit):
        dx = np.linalg.inv(Q).dot(b -A.dot(x))
        x = x + lamb*dx
        if np.linalg.norm(dx)<tol:
            break
    print("x:\n", x) 
    return x

gseidel(A,b)<jupyter_output>x:
 [ 2.  1.]<jupyter_text>Here, we have incorporated a so-called over-relaxation parameter, lambda. Instead of using $x + dx$, we use $x + lamb*x$ to compute the next iterate.a value of lambda between 1 and 2 will accelerate convergence of the Gauss-Seidel algorithm.  gauss seidel methodhttp://austingwalters.com/gauss-seidel-method/Using the Gauss-Seidel MethodThe method is fairly straight forward, given a standard system of linear equations, Ax = b. Where, A is a matrix (often representing a series of equations), x is a vector of x variables (Gauss-Seidel method is used to solve this vector) and b is the solution vector. In Gauss-Seidel method, we then split the A matrix into Upper (U) and Lower (L) matrices (the lower matrix in this case also contains the diagonal), then iterate using the following method:![](http://108.61.119.12/wp-content/uploads/2014/05/gauss-equation.png)<jupyter_code>"""
author:  hahnicity
https://github.com/hahnicity/ace/blob/master/ace/solvers.py
Solve a simple linear equation Ax = b in a particular way
"""


def gauss_seidel(A, b, maxit=1000, x0=None):
    """
    Solve a linear equation by the gauss seidel iteration outlined in the book
    Follows the eq:
        x = inv(L)*(b - U*x)
    """
    l = tril(A)
    upper_plus_diagonal = triu(A)
    u = upper_plus_diagonal - np.diag(np.diag(A))
    # If we have not provided an initial array for x make a new one
    if x0 == None:
        x = array([1 for _ in range(A.shape[1])])
    else:
        x = x0

    for _ in np.arange(maxit):
        x = inv(l).dot(b - u.dot(x))
    return x




gauss_seidel(A, b)<jupyter_output><empty_output><jupyter_text>The following script solves a linear equation using Gauss-Seideliteration with default value of 1 for the over-relaxation parameter:<jupyter_code>A = np.array([[3, 1],
               [2, 5]])
b = np.array([7 , 9])
x = gauss_seidel(A,b,maxit)
print("x:\n",x)<jupyter_output>x:
 [ 2.  1.]<jupyter_text>When A=[3 2; 4 1], however,the algorithm diverges.<jupyter_code>A = np.array([[3, 2],
               [4, 1]])
b = np.array([7 , 9])
x = gauss_seidel(A,b,maxit)
print("x:\n",x)<jupyter_output>x:
 [ nan  nan]
{'fraction_non_alphanumeric': 0.11462888820449646, 'fraction_numerical': 0.03578688019710502, 'mean_word_length': 3.6295979469632162, 'pattern_counts': {'":': 0, '<': 67, '<?xml version=': 0, '>': 64, 'https://': 4, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '279356', 'license': 'MIT', 'n_tokens_mistral': 5781, 'n_tokens_neox': 5313, 'n_words': 1989, 'path': 'Chapter02.ipynb', 'repo_name': 'lnsongxf/Applied_Computational_Economics_and_Finance', 'size': '16235'}