Algorithms Homework 1#

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

Linear Algebra Review#

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.

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.

# Add your solution here

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

b = 
 [6 8 4]

# Add 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]]

# Add your solution here

[1. 1. 1.]

Solve using LU decomposition#

Do the following:

  1. Use linalg.lu(A) to calculate \(P\), \(L\), and \(U\).

  2. 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.

  3. Use your function to solve the linear system.

# 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.

# 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 = []
    
    # Add your solution here
    
    return x

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

Solve using linalg.solve#

# Add your solution here

x =  [1. 1. 1.]

Eigenvalues#

Calculate eigenvalues by hand (on paper)#

Calculate the eigenvalues for the following matrix:

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

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.

Calculate eigenvalues using linalg.eig#

Calculate the eigenvalues and corresponding (right) eigenvectors.

# Add your solution here

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

# Add 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]]

Definiteness#

Based on your calculations above, is this matrix

  1. positive definite

  2. positive semi definite

  3. negative definite

  4. negative semi definite

  5. indefinite or

  6. 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.

Singular Value Decomposition#

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

SVD calculation using linalg.svd#

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

\[\begin{split} A_3 = \begin{bmatrix} 1 & 2 & 0 & -1 \\ 2 & 4 & 10^{-6} & -2 \\ 0 & 2 & -1 & 10^{-3} \end{bmatrix} \end{split}\]
# Add 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]]

# Add 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]]

Condition number#

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

  1. Using SVD results

  2. Using np.linalg.cond()

# Add your solution here

Condition number =  16955285.123025477

# Add your solution here

Condition number =  16955285.123025477

Linear system#

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

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

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

# Add your solution here

Condition number =  181.90540228772272

Answer:

Make it singular#

Do/answer the following:

  1. Propose a change to a single entry in \(A_4\) to make it singular.

  2. What are the singular values and condition number after this change?

  3. What can you say about the solution to \(A_4 \cdot x=b\) after the change?

Answer:

Convexity#

Please turn in via Gradescope.

Determine if the following functions are convex#

  1. \(x^2 + ax + b, ~x \in \mathbb{R}\)

  2. \(x^3, ~x \in \mathbb{R}\)

  3. \(| x |, ~ x \in [ -5, 5]\)

  4. \(x + y, ~ x \in \mathbb{R}, ~ y \in \mathbb{R}\)

  5. \(x \cdot y, ~ x \in \mathbb{R}, ~ y \in \mathbb{R}\)

  6. \(\log(x), ~ x \in (0,1]\)

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\)

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.

Convexity implies PSD#

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

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.