.. _SVDHW:

Singular Value Decomposition Homework
=====================================

#. Find the singular values of the following matrix and determine the
   matrix’s rank from the singular values. Show the Python commands.

   .. math::

      \mathbf{A} = \begin{bmatrix}
                      1  & -2 & -4 & 8 \\ 
                      14  & 4 & 12 & 12 \\ 
                      7  & 8 & 1 & 5 \\ 
                      -2  & 1 & -8 & 3
                   \end{bmatrix}

   ::

    A = np.array([[1, -2, -4, 8], [14, 4, 12, 2], [7, 8, 1, 5], [-2, 1, -8, 3]]).astype(float)



#. Find the singular values of the following matrix. Find the
   condition number of :math:`\bf{A}` from the ratio of the largest
   to smallest singular values (:math:`\sigma_1/\sigma_3`). Verify
   your answer with the ``np.linalg.cond`` function.

   .. math::

      \mathbf{A} = \begin{bmatrix}
                      -8  & 13 & 3 \\ 
                      -1  & 11 & 5 \\ 
                      -7  & 14 & 4
                   \end{bmatrix}

   ::

    A = np.array([[-8, 13, 3], [-1, 11, 5], [-7, 14, 4]]).astype(float)

#. Write a function to compute the rank of a matrix from its singular
   values. Set a threshold value of :math:`10^{-13}` (1e-13) for 
   determining if the singular
   values should be counted as zero or nonzero values.

   *Hint*: Create a logical array and count the nonzero values with
   ``np.count_nonzero``.

   Use the following code as a starting point. Fill in the code for the
   function. The testing code is provided.

   ::

    # -*- coding: utf-8 -*-
    """
    myrank

    Find the rank of a matrix from the singular values
    """

    import numpy as np
    from scipy.linalg import svd

    def myrank(A: np.array) -> int:

    if __name__ == '__main__':
        A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]).astype(float)
        U,s,Vt = svd(A)
        S = np.diag(s)
        S[2,2] = 1e-14
        print( myrank(U@S@Vt) == 2)
        S[2,2] = 1e-12
        print( myrank(U@S@Vt) == 3)
        
   
#. Write functions to replicate Python's ``scipy.linalg.orth`` function from
   the SVD. The function should return the orthogonal basis for a
   matrix.  The number of columns returned should be determined by the
   matrix's rank.

   Since the basis matrix comes from the SVD, it will be orthogonal, so we don't
   need to test for that. We should test that the correct number of columns
   is returned and that the columns form a basis for :math:`\bf{A}`, 
   which means that :math:`\bf{A}` is in the column space for the 
   basis (:math:`\bf{U}`). There are two ways to test that
   :math:`\bf{A}` is in column space of  :math:`\bf{U}`. 
   The rank augmented matrix :math:`[\mathbf{U A}]` should be the same
   as the rank of :math:`\bf{U}`. The other test is the projection of
   :math:`\bf{A}` onto :math:`\bf{U}` should yield exactly :math:`\bf{A}`.
   The testing code below uses the second test.

   Use the following code as a starting point. Fill in the code for the
   function. The testing code is provided.

   ::

    # -*- coding: utf-8 -*-
    """
    myorth
    Find the orthogonal basis for Col(A) from the SVD
    """

    import numpy as np
    from scipy.linalg import svd

    def myorth(A: np.array) -> np.array:

    if __name__ == '__main__':
        # square, singular matrix
        A = np.random.randint(-5, 30, (4,4)).astype(float)
        A[:,[2]] = A[:,[0]]/2 + A[:,[1]]/2
        A[:,[3]] = A[:,[0]]/3 + A[:,[1]]/3
        U = myorth(A)
        # Col(A) should have 2 columns
        m,n = U.shape
        print(n == 2)
        # Is A in column space of U?
        # The projection of A onto U should yield A
        # if U is a basis for Col(A)
        projection = U @ (U.T @ A)
        print(np.allclose(projection, A))
        # over-determined
        A = np.random.randint(-5, 30, (5,2)).astype(float)
        U = myorth(A)
        # Is A in column space of U?
        projection = U @ (U.T @ A)
        print(np.allclose(projection, A))
        # Under-determined
        A = np.random.randint(-5, 30, (3,5)).astype(float)
        U = myorth(A)
        # Is A in column space of U?
        projection = U @ (U.T @ A)
        print(np.allclose(projection, A))

#. Write functions to replicate Python's 
   ``scipy.linalg.null_space``
   function from the SVD. The function should return
   the null space of the matrix. The number of columns returned by the
   function should be determined by the matrix's rank.

   Use the following code as a starting point. Fill in the code for the
   function. The testing code is provided.

   ::

    # -*- coding: utf-8 -*-
    """
    mynull
    Find the null space of A (v such that Av = 0). A must be singular.
    """

    import numpy as np
    from scipy.linalg import svd

    def mynull(A: np.array) -> np.array:

    if __name__ == '__main__':
        # square matrix
        A = np.random.randint(-5, 30, (4,4)).astype(float)
        A[:,[2]] = A[:,[0]] + A[:,[1]]/2
        V = mynull(A)
        print(np.allclose(A@V, np.zeros((4,1))))
        # print(V)