This notebook contains material from CBE60499; content is available on Github.

< 3.6 Descent and Globalization | Contents | Tag Index | 3.8 Algorithms Homework 2 >

3.7 Algorithms Homework 1¶

In [1]:

# Load required Python libraries.
import matplotlib.pyplot as plt
import numpy as np
from scipy import linalg

3.7.1 Linear Algebra Review¶

3.7.1.1 Exact solution¶

Solve the following system of linear equations BY HAND with exact arithmetic:

$$3 x_1 + 2 x_2 + x_3 = 6$$$$-x_1 + 4 x_2 + 5x_3 = 8$$$$2x_1 -8 x_2 + 10 x_3 = 4$$

Turn in your work via Gradescope.

3.7.1.2 Solve using the inverse¶

Solve the same linear system by first inverting the matrix A and performing matrix multiplication. You should use Python and SciPy.

Tutorials and documentation: https://docs.scipy.org/doc/scipy-1.1.0/reference/tutorial/linalg.html

In [2]:

# YOUR SOLUTION HERE

A = 
 [[ 3  2  1]
 [-1  4  5]
 [ 2 -8 10]]

b = 
 [6 8 4]

In [3]:

# YOUR SOLUTION HERE

[[ 2.85714286e-01 -1.00000000e-01  2.14285714e-02]
 [ 7.14285714e-02  1.00000000e-01 -5.71428571e-02]
 [ 4.62592927e-19  1.00000000e-01  5.00000000e-02]]

In [4]:

# YOUR SOLUTION HERE

[1. 1. 1.]

3.7.1.3 Solve using LU decomposition¶

Do the following:

Use linalg.lu(A) to calculate $P$, $L$, and $U$.
Write a function to solve any linear system given the factorization and $b$. You may not use linalg.solve in your function. Instead, you should write loops to implement back substitution.
Use your function to solve the linear system.

In [5]:

# You need to define the matrix A1 somewhere (or change the variable name)
P, L, U = linalg.lu(A1)

# Permutation matrix
print(P)

# Lower diagonal matrix
print(L)

# Upper diagonal matrix
print(U)

[[1. 0. 0.]
 [0. 0. 1.]
 [0. 1. 0.]]
[[ 1.          0.          0.        ]
 [ 0.66666667  1.          0.        ]
 [-0.33333333 -0.5         1.        ]]
[[ 3.          2.          1.        ]
 [ 0.         -9.33333333  9.33333333]
 [ 0.          0.         10.        ]]

Tip: We discussed this algorithm in class. First write pseudocode on paper (how to translate our notes into loops, logical statements, etc.?). Once you are happy with the logic, code it up.

In [6]:

# Write function here.
def my_lu_solve(P, L, U, b, LOUD=True):
    ''' 
    Solves linear system Ax = b given PLU decomposition of A
    
    Arguments:
        P - N by N permutation matrix
        L - N by N lower triangular matrix with 1 on diagonal
        U - N by N upper triangular matrix
        b - N by 1 vector
    
    Returns:
        x - N by 1 solution to linear system
    '''
    
    # Define x so this does not give an error. You can delete the line below.
    x = []
    
    # YOUR SOLUTION HERE
    
    return x

In [7]:

x = my_lu_solve(P, L, U, b, LOUD=True)

Forward solve...
y[ 0 ]= 6.0
y[ 1 ]= 0.0
y[ 2 ]= 10.0

Backward solve...
x[ 2 ]= 1.0
x[ 1 ]= 1.0
x[ 0 ]= 1.0

3.7.1.4 Solve using `linalg.solve`¶

In [8]:

# YOUR SOLUTION HERE

x =  [1. 1. 1.]

3.7.2 Eigenvalues¶

3.7.2.1 Calculate eigenvalues by hand (on paper)¶

Calculate the eigenvalues for the following matrix:

$$ A = \begin{bmatrix} 0 & -4 & -6 \\ -1 & 0 & -3 \\ 1 & 2 & 5 \end{bmatrix} $$

You may need to do this for a small system on an exam to characterize stationary points of an optimization problem. Hence I would like you to practice at least once on the homework.

3.7.2.2 Calculate eigenvalues using `linalg.eig`¶

Calculate the eigenvalues and corresponding (right) eigenvectors.

In [9]:

# YOUR SOLUTION HERE

A2 = 
 [[ 0 -4 -6]
 [-1  0 -3]
 [ 1  2  5]]

In [10]:

# YOUR SOLUTION HERE

Eigenvalues =  [1.+0.00000000e+00j 2.+1.63168795e-15j 2.-1.63168795e-15j]
Eigenvectors = 
[[-0.81649658+0.j         -0.77616123+0.j         -0.77616123-0.j        ]
 [-0.40824829+0.j         -0.26081756-0.31398875j -0.26081756+0.31398875j]
 [ 0.40824829+0.j          0.43259879+0.20932583j  0.43259879-0.20932583j]]

3.7.2.3 Definiteness¶

Based on your calculations above, is this matrix

positive definite
positive semi definite
negative definite
negative semi definite
indefinite or
cannot say without additional calculations?

Briefly comment to justify your answer.

Note: In this class, "briefly comment" means write a few sentences, sketch a picture, write an equation or some combination... whichever you feel is necessary to succinctly provide justification.

3.7.2.4 Singular Value Decomposition¶

TODO: Make this a separate problem (after the assignment is turned in this year).

3.7.2.5 SVD calculation using `linalg.svd`¶

Calculate the singular value decomposition of the following matrix using linalg.svd:

$$ A_3 = \begin{bmatrix} 1 & 2 & 0 & -1 \\ 2 & 4 & 10^{-6} & -2 \\ 0 & 2 & -1 & 10^{-3} \end{bmatrix} $$

In [11]:

# YOUR SOLUTION HERE

A = 
 [[ 1.e+00  2.e+00  0.e+00 -1.e+00]
 [ 2.e+00  4.e+00  1.e-06 -2.e+00]
 [ 0.e+00  2.e+00 -1.e+00  1.e-03]]

In [12]:

# YOUR SOLUTION HERE

U = 
 [[-4.25830830e-01 -1.36631105e-01 -8.94427217e-01]
 [-8.51661650e-01 -2.73262656e-01  4.47213544e-01]
 [-3.05516837e-01  9.52186674e-01  1.91553449e-07]]
S = 
 [5.73316009e+00 1.45975216e+00 3.38134101e-07]
Vh = 
 [[-3.71375314e-01 -8.49329490e-01  5.32892822e-02  3.71322025e-01]
 [-4.67994798e-01  3.68597169e-01 -6.52293569e-01  4.68647091e-01]
 [-3.77573248e-01  3.77856348e-01  7.56090836e-01  3.78139751e-01]
 [ 7.07460026e-01 -3.53377112e-04 -9.52352187e-10  7.06753272e-01]]

3.7.2.6 Condition number¶

What is the condition number of the matrix? Calculate it two ways:

Using SVD results
Using np.linalg.cond()

In [13]:

# YOUR SOLUTION HERE

Condition number =  16955285.123025477

In [14]:

# YOUR SOLUTION HERE

Condition number =  16955285.123025477

3.7.2.7 Linear system¶

Consider the linear system $A_4 \cdot x = b$ with $b = [1, 0, 3]^T$ and

$$ A_4 = \begin{bmatrix} 1 & 2 & 0 \\ 2 & 4 & 10^{-1} \\ 0 & 2 & -1 \end{bmatrix} ~. $$

Approximately how much uncertainty $||\Delta b||_2$ can you tolerate if you want the uncertainty $||\Delta x||_2$ to be less than $10^{-2}$?

In [15]:

# YOUR SOLUTION HERE

Condition number =  181.90540228772272

Answer:

3.7.2.8 Make it singular¶

Do/answer the following:

Propose a change to a single entry in $A_4$ to make it singular.
What are the singular values and condition number after this change?
What can you say about the solution to $A_4 \cdot x=b$ after the change?

Answer:

3.7.3 Convexity¶

Please turn in via Gradescope.

3.7.3.1 Determine if the following functions are convex¶

$x^2 + ax + b, ~x \in \mathbb{R}$
$x^3, ~x \in \mathbb{R}$
$| x |, ~ x \in [ -5, 5]$
$x + y, ~ x \in \mathbb{R}, ~ y \in \mathbb{R}$
$x \cdot y, ~ x \in \mathbb{R}, ~ y \in \mathbb{R}$
$\log(x), ~ x \in (0,1]$

3.7.3.2 Prove the following properties¶

Consider twice differentiable function $f(x): \mathbb{R}^n \rightarrow \mathbb{R}$.

Recall that $f(x)$ is convex on $x \in \mathbb{R}^n$ if and only if

$f(x+p) \geq f(x) + \nabla f(x)^T p$ for all $x, p \in \mathbb{R}^n$

3.7.3.2.1 PSD implies Convexity¶

Prove that if $\nabla^2 f(x)$ is positive semidefinite for all $x \in \mathbb{R}^n$, then $f(x)$ must be convex.

3.7.3.2.2 Convexity implies PSD¶

Prove that if $f(x)$ is convex then $\nabla^2 f(x)$ must be positive semidefinite.

3.7.3.2.3 PD implies Strictly Convex¶

Prove that if $\nabla^2 f(x)$ is positive definite for all $x \in \mathbb{R}^n$, then $f(x)$ must be strictly convex.