#! /usr/bin/env python3
#
def svd_product ( ):

## svd_product demonstrates that the SVD factors reproduce the matrix.
#
  import numpy as np

  print ( '' )
  print ( 'svd_product' )
  print ( '  Show that the product of the SVD factors is the original matrix.' )
#
#  A is a square matrix.
#
  m = 5
  n = 5
  A = np.array ( [
    [17, 24,  1,  8, 15],
    [23,  5,  7, 14, 16],
    [ 4,  6, 13, 20, 22],
    [10, 12, 19, 21,  3],
    [11, 18, 25,  2,  9] ] )
  print ( '' )
  print ( 'A:')
  print ( A )
  U, svec, V = np.linalg.svd ( A )
  S = np.diag ( svec )               # This only works when A is square.
  SV = np.matmul ( S, V )
  USV = np.matmul ( U, SV )
  diff = np.linalg.norm ( A - USV )
  print ( '' )
  print ( '  || A - U*S*V|| = ', diff )
#
# B is a "tall" matrix
#
  m = 6
  n = 4
  B = np.array ( [
    [17, 24,  1,  8],
    [23,  5,  7, 14],
    [ 4,  6, 13, 20],
    [10, 12, 19, 21],
    [11, 18, 25,  2],
    [15, 16, 22,  3],
    [ 2,  0,  1,  9] ] )
  print ( '' )
  print ( 'B:')
  print ( B )
  U, svec, V = np.linalg.svd ( B )
  S = np.zeros ( B.shape )
  for i in range ( 0, min ( m, n ) ):
    S[i,i]= svec[i]
  SV = np.matmul ( S, V )
  USV = np.matmul ( U, SV )
  diff = np.linalg.norm ( B - USV )
  print ( '' )
  print ( '  || B - U*S*V|| = ', diff )
#
# C is a "wide" matrix
#
  m = 4
  n = 6
  C = np.array ( [
    [17, 24,  1,  8, 23,  5],
    [ 7, 14,  4,  6, 13, 20],
    [10, 12, 19, 21, 11, 18],
    [25,  2, 15, 16, 22,  3] ] )
  print ( '' )
  print ( 'C:')
  print ( C )
  U, svec, V = np.linalg.svd ( C )
  S = np.zeros ( C.shape )
  for i in range ( 0, min ( m, n ) ):
    S[i,i]= svec[i]
  SV = np.matmul ( S, V )
  USV = np.matmul ( U, SV )
  diff = np.linalg.norm ( C - USV )
  print ( '' )
  print ( '  || C - U*S*V|| = ', diff )

  return

if ( __name__ == '__main__' ):
  svd_product ( )