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

  import numpy as np

  print ( '' )
  print ( 'svd_approx' )
  print ( '  Show that A can be approximated by partial SVD products.' )
#
#  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 )
#
#  Compute A0, the column averages, and subtract it.
#  A = A0 + Ar
#
  u0 = np.ones ( m )
  v0 = np.mean ( A, axis = 0 )
  A0 = np.outer ( u0, v0 )
#
#  Get the singular value decomposition of the remainder.
#
  U, s, V = np.linalg.svd ( A - A0 )

  diff = np.linalg.norm ( A )
  print ( '' )
  print ( '  Norm of A = %g' % ( diff ) )
#
#  Compute a 0 component approximation to A.
#
  diff = np.linalg.norm ( A - A0 )
  print ( '' )
  print ( '  Norm of A-A0 = %g' % ( diff ) )
  print ( '  Square root of sum of singular values = ', np.sqrt ( np.sum ( s**2 ) ) )
#
#  Compute a 1 component approximation to A.
#
  A1 = s[0] * np.outer ( U[:,0], V[0,:] )
  diff = np.linalg.norm ( A - A0 - A1 )
  print ( '' )
  print ( '  Norm of A-A0-A1 = %g' % ( diff ) )
  print ( '  Singular value percentage = ', \
    100 * np.sqrt ( np.sum ( s[0:1]**2 ) ) / np.sqrt ( np.sum ( s**2 ) ) )
  print ( '' )
  print ( 'A0+A1:' )
  print ( A0+A1 )
#
#  Compute a 2 component approximation to A.
#
  A2 = s[1] * np.outer ( U[:,1], V[1,:] )
  diff = np.linalg.norm ( A - A0 - A1 - A2 )
  print ( '' )
  print ( '  Norm of A-A0-A1-A2 = %g' % ( diff ) )
  print ( '  Singular value percentage = ', \
    100 * np.sqrt ( np.sum ( s[0:2]**2 ) ) / np.sqrt ( np.sum ( s**2 ) ) )
  print ( '' )
  print ( 'A0+A1+A2:' )
  print ( A0+A1+A2 )
#
#  Compute a 3 component approximation to A.
#
  A3 = s[2] * np.outer ( U[:,2], V[2,:] )
  diff = np.linalg.norm ( A - A0 - A1 - A2 - A3 )
  print ( '' )
  print ( '  Norm of A-A0-A1-A2-A3 = %g' % ( diff ) )
  print ( '  Singular value percentage = ', \
    100 * np.sqrt ( np.sum ( s[0:3]**2 ) ) / np.sqrt ( np.sum ( s**2 ) ) )
#
#  Compute a 4 component approximation to A.
#
  A4 = s[3] * np.outer ( U[:,3], V[3,:] )
  diff = np.linalg.norm ( A - A0 - A1 - A2 - A3 - A4 )
  print ( '' )
  print ( '  Norm of A-A0-A1-A2-A3-A4 = %g' % ( diff ) )
  print ( '  Singular value percentage = ', \
    100 * np.sqrt ( np.sum ( s[0:4]**2 ) ) / np.sqrt ( np.sum ( s**2 ) ) )

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