function svd_test ( )

%*****************************************************************************80
%
%% svd_test() tests svd().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    17 June 2021
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'svd_test\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  svd() computes the singular value decomposition.\n' );

  svd_demo ( 5, 3 );
  svd_demo ( 5, 5 );
  svd_demo ( 5, 8 );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'svd_test():\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function svd_demo ( m, n )

%*****************************************************************************80
%
%% svd_demo() demonstrates the singular value decomposition.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modififed:
%
%    27 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Usage:
%
%    svd_test ( m, n )
%
%  Input:
%
%    integer M, N, the number of rows and columns
%    of the matrix.  If M and N are not provided on the user
%    is prompted to supply them.
%
%  Local:
%
%    real A(M,N), the matrix whose singular value
%    decomposition we are investigating.
%
%    real S(M,N), the diagonal factor
%    int the singular value decomposition of A.
%
%    real U(M,M), the first orthogonal factor
%    in the singular value decomposition of A.
%
%    real V(N,N), the second orthogonal factor
%    in the singular value decomposition of A.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'svd_demo():\n' );
  fprintf ( 1, '  Test the SVD for a %d by %d matrix\n', n, n );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  A real MxN matrix A can be factored as:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    A = U * S * V''\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  where\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    U = MxM orthogonal,\n' );
  fprintf ( 1, '    S = MxN zero except for diagonal,\n' );
  fprintf ( 1, '    V = NxN orthogonal.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The diagonal of S contains only nonnegative numbers\n' );
  fprintf ( 1, '  and these are arranged in descending order.\n' );
%
%  If at least one command line argument, it's M.
%
  if ( nargin < 1 )
    m = input ( 'Enter the value of M (M <= 10 is best):' );
  end
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix row order    M = %d\n', m );
%
%  If at least two command line argument, it's N.
%
  if ( nargin < 2 )
    n = input ( 'Enter the value of N (N <= 10 is best):' );
  end
  fprintf ( 1, '  Matrix column order N = %d\n', n );
%
%  Generate the matrix A.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Choose a "random" matrix A, with integral values.\n' );

  a = round ( 10.0 * randn ( m, n ) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The matrix A:\n' );
  disp ( a );
%
%  Get the SVD, as MATLAB computes it.
%
  [ u, s, v ] = svd ( a );
%
%  Print the SVD.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Orthogonal factor U:\n' );
  disp ( u );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Diagonal factor S:\n' );
  disp ( s );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Orthogonal factor V:\n' );
  disp ( v );
%
%  Check that A = U * S * V'.
%
  svd_product_test ( m, n, a, u, s, v );
%
%  Compute the norm of the difference between A and the successive
%  sums of rank one approximants.
%
  rank_one_test ( m, n, a, u, s, v );
%
%  Actually print the sums of rank one approximants.
%
  rank_one_print_test ( m, n, a, u, s, v );
%
%  Compute the pseudoinverse.
%
  a_pseudo = pseudo_inverse ( m, n, u, s, v );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Pseudoinverse of A:\n' );
  disp ( a_pseudo );
%
%  Test A*A+ = I+, A+*A = I+
%
  pseudo_product_test ( m, n, a, a_pseudo );
%
%  Demonstrate the use of the pseudoinverse for linear systems.
%
  pseudo_linear_solve_test ( m, n, a, a_pseudo );
  
  return
end
function a_pseudo = pseudo_inverse ( m, n, u, s, v )

%*****************************************************************************80
%
%% pseudo_inverse() computes the pseudoinverse.
%
%  Discussion:
%
%    Given the singular value decomposition of a real MxN matrix A:
%
%      A = U * S * V'
%
%    where 
%
%      U is MxM orthogonal,
%      S is MxN, and entirely zero except for the diagonal;
%      V is NxN orthogonal.
%
%    the pseudo inverse is the NxM matrix A+ with the form
%
%      A+ = V * S+ * U'
%
%    where 
%
%      S+ is the NxM matrix whose nonzero diagonal elements are
%      the inverses of the corresponding diagonal elements of S.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    13 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer M, N, the number of rows and columns in the matrix A.
%
%    real A(M,N), the matrix whose singular value
%    decomposition we are investigating.
%
%    real U(M,M), S(M,N), V(N,N), the factors
%    that form the singular value decomposition of A.
%
%  Output:
%
%    real A_PSEUDO(N,M), the pseudoinverse of A.
%
  sp(1:n,1:m) = 0.0;
  for i = 1 : min ( m, n )
    if ( s(i,i) ~= 0.0 )
      sp(i,i) = 1.0 / s(i,i);
    end
  end

  a_pseudo(1:n,1:m) = v(1:n,1:n) * sp(1:n,1:m) * u(1:m,1:m)';

  return
end
function pseudo_linear_solve_test ( m, n, a, a_pseudo )

%*****************************************************************************80
%
%% pseudo_linear_solve_test() uses the pseudoinverse for linear systems.
%
%  Discussion:
%
%    Given an MxN matrix A, and its pseudoinverse A+:
%
%      "Solve" A  * x = b by x = A+  * b.
%
%      "Solve" A' * x = b by x = A+' * b.
%
%    When the system is overdetermined, the solution minimizes the
%    L2 norm of the residual.  
%
%    When the system is underdetermined, the solution
%    is the solution of minimum L2 norm.     
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    14 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer M, N, the number of rows and columns in the matrix A.
%
%    real A(M,N), the matrix whose singular value
%    decomposition we are investigating.
%
%    real A_PSEUDO(N,M), the pseudo_inverse of A.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'pseudo_linear_solve_test():\n' );
%
%  A * x = b, b in range of A.
%
  xn1(1:n,1) = round ( 10.0 * randn ( n, 1 ) );
  bm(1:m,1) = a(1:m,1:n) * xn1(1:n,1);
  xn2(1:n,1) = a_pseudo(1:n,1:m) * bm(1:m,1);
  rm(1:m,1) = a(1:m,1:n) * xn2(1:n,1) - bm(1:m,1);

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Given:\n' );
  fprintf ( 1, '    b = A * x1\n' );
  fprintf ( 1, '  so that b is in the range of A, solve\n' );
  fprintf ( 1, '    A * x = b\n' );
  fprintf ( 1, '  using the pseudoinverse:\n' );
  fprintf ( 1, '    x2 = A+ * b.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Norm of x1 = %f\n', norm ( xn1 ) );
  fprintf ( 1, '  Norm of x2 = %f\n', norm ( xn2 ) );
  fprintf ( 1, '  Norm of residual = %f\n', norm ( rm ) );
%
%  A * x = b, b not in range of A.
%
  if ( n < m )

    fprintf ( 1, '\n' );
    fprintf ( 1, '  For N < M, most systems A*x=b will not be\n' );
    fprintf ( 1, '  exactly and uniquely solvable, except in the\n' );
    fprintf ( 1, '  least squares sense.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Here is an example:\n' );

    b8(1:m,1) = randn ( m, 1 );
    xn2(1:n,1) = a_pseudo(1:n,1:m) * bm(1:m,1);
    rm(1:m,1) = a(1:m,1:n) * xn2(1:n,1) - bm(1:m,1);

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Given b is NOT in the range of A, solve\n' );
    fprintf ( 1, '    A * x = b\n' );
    fprintf ( 1, '  using the pseudoinverse:\n' );
    fprintf ( 1, '    x2 = A+ * b.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Norm of x2 = %f\n', norm ( xn2 ) );
    fprintf ( 1, '  Norm of residual = %f\n', norm ( rm ) );

  end
%
%  A' * x = b, b is in the range of A'.
%
  xm1(1:m,1) = round ( 10.0 * randn ( m, 1 ) );
  xm1(1:m,1) = round ( 10.0 * xm1(1:m,1) );
  bn(1:n,1) = a(1:m,1:n)' * xm1(1:m,1);
  xm2(1:m,1) = a_pseudo(1:n,1:m)' * bn(1:n,1);
  rn(1:n,1) = a(1:m,1:n)' * xm2(1:m,1) - bn(1:n,1);

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Given:\n' );
  fprintf ( 1, '    b = A'' * x1\n' );
  fprintf ( 1, '  so that b is in the range of A'', solve\n' );
  fprintf ( 1, '    A'' * x = b\n' );
  fprintf ( 1, '  using the pseudoinverse:\n' );
  fprintf ( 1, '    x2 = A+'' * b.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Norm of x1 = %f\n', norm ( xm1 ) );
  fprintf ( 1, '  Norm of x2 = %f\n', norm ( xm2 ) );
  fprintf ( 1, '  Norm of residual = %f\n', norm ( rn ) );
%
%  A' * x = b, b is not in the range of A'.
%
  if ( m < n )

    fprintf ( 1, '\n' );
    fprintf ( 1, '  For M < N, most systems A''*x=b will not be\n' );
    fprintf ( 1, '  exactly and uniquely solvable, except in the\n' );
    fprintf ( 1, '  least squares sense.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Here is an example:\n' );

    bn(1:n,1) = randn ( n, 1 );
    xm2(1:m,1) = a_pseudo(1:n,1:m)' * bn(1:n,1);
    rn(1:n,1) = a(1:m,1:n)' * xm2(1:m,1) - bn(1:n,1);

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Given b is NOT in the range of A'', solve\n' );
    fprintf ( 1, '    A'' * x1 = b\n' );
    fprintf ( 1, '  using the pseudoinverse:\n' );
    fprintf ( 1, '    x2 = A+'' * b.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Norm of x2 = %f\n', norm ( xm2 ) );
    fprintf ( 1, '  Norm of residual = %f\n', norm ( rn ) );

  end

  return
end
function pseudo_product_test ( m, n, a, a_pseudo )

%*****************************************************************************80
%
%% pseudo_product_test() examines the products A*A+ and A+*A.
%
%  Discussion:
%
%    Given an MxN matrix A, and its pseudoinverse A+, we must have
%
%      A+ * A * A+ = A+
%      A * A+ * A = A
%      ( A * A+ )' = A * A+ (MxM symmetry)
%      ( A+ * A )' = A+ * A (NxN symmetry)
%
%    If M <= N, A * A+ may be "interesting" (equal to or "like" the identity),
%    if N <= M, A+ * A may be "interesting" (equal to or "like" the identity).
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    14 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer M, N, the number of rows and columns in the matrix A.
%
%    real A(M,N), the matrix whose singular value
%    decomposition we are investigating.
%
%    real A_PSEUDO(N,M), the pseudo_inverse of A.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'pseudo_product_test()\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The following relations MUST hold:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   A  * A+ * A  = A\n' );
  fprintf ( 1, '   A+ * A  * A+ = A+\n' );
  fprintf ( 1, ' ( A  * A+ ) is MxM symmetric;\n' );
  fprintf ( 1, ' ( A+ * A  ) is NxN symmetric\n' );

  bmn(1:m,1:n) = a(1:m,1:n) * a_pseudo(1:n,1:m) * a(1:m,1:n);

  dif1 = norm ( a - bmn, 'fro' );

  bnm(1:n,1:m) = a_pseudo(1:n,1:m) * a(1:m,1:n) * a_pseudo(1:n,1:m);

  dif2 = norm ( a_pseudo - bnm, 'fro' );

  bmm(1:m,1:m) = a(1:m,1:n) * a_pseudo(1:n,1:m);

  dif3 = norm ( bmm - transpose ( bmm ), 'fro' );

  bnn(1:n,1:n) = a_pseudo(1:n,1:m) * a(1:m,1:n);

  dif4 = norm ( bnn - transpose ( bnn ), 'fro' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Here are the Frobenius norms of the errors\n' );
  fprintf ( 1, '  in these relationships:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   A  * A+ * A  = A            %f\n', dif1 );
  fprintf ( 1, '   A+ * A  * A+ = A+           %f\n', dif2 );
  fprintf ( 1, ' ( A  * A+ ) is MxM symmetric; %f\n', dif3 );
  fprintf ( 1, ' ( A+ * A  ) is NxN symmetric; %f\n', dif4 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  In some cases, the matrix A * A+\n' );
  fprintf ( 1, '  may be interesting (if M <= N, then\n' );
  fprintf ( 1, '  it MIGHT look like the identity.)\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  A * A+:\n' );
  disp ( bmm );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  In some cases, the matrix A+ * A\n' );
  fprintf ( 1, '  may be interesting (if N <= M, then\n' );
  fprintf ( 1, '  it MIGHT look like the identity.)\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  A+ * A:\n' );
  disp ( bnn );

  return
end
function rank_one_print_test ( m, n, a, u, s, v )

%*****************************************************************************80
%
%% rank_one_print_test() prints the sums of rank one matrices.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    13 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer M, N, the number of rows and columns in the matrix A.
%
%    real A(M,N), the matrix whose singular value
%    decomposition we are investigating.
%
%    real U(M,M), S(M,N), V(N,N), the factors
%    that form the singular value decomposition of A.
%
  a_norm = norm ( a, 'fro' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'rank_one_print_test():\n' );
  fprintf ( 1, '  Print the sums of R rank one matrices.\n' );

  for r = 0 : min ( m, n )

    usv(1:m,1:n) = u(1:m,1:r) * s(1:r,1:r) * v(1:n,1:r)';

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Rank R = %d:\n', r );
    disp ( usv );

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Original matrix A:\n' );
  disp ( a );

  return
end
function rank_one_test ( m, n, a, u, s, v )

%*****************************************************************************80
%
%% rank_one_test() compares A to the sum of rank one matrices.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    13 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer M, N, the number of rows and columns in the matrix A.
%
%    real A(M,N), the matrix whose singular value
%    decomposition we are investigating.
%
%    real U(M,M), S(M,N), V(N,N), the factors
%    that form the singular value decomposition of A.
%
  a_norm = norm ( a, 'fro' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'rank_one_test():\n' );
  fprintf ( 1, '  Compare A to the sum of R rank one matrices.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '         R    Absolute      Relative\n' );
  fprintf ( 1, '              Error         Error\n' );
  fprintf ( 1, '\n' );

  for r = 0 : min ( m, n )

    usv(1:m,1:n) = u(1:m,1:r) * s(1:r,1:r) * v(1:n,1:r)';

    dif_norm = norm ( a - usv, 'fro' );

    fprintf ( 1, '  %8d  %14f  %14f\n', r, dif_norm, dif_norm / a_norm );

  end

  return
end
function svd_product_test ( m, n, a, u, s, v )

%*****************************************************************************80
%
%% svd_product_test() tests that A = U * S * V'.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    13 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer M, N, the number of rows and columns in the matrix A.
%
%    real A(M,N), the matrix whose singular value
%    decomposition we are investigating.
%
%    real U(M,M), S(M,N), V(N,N), the factors
%    that form the singular value decomposition of A.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'svd_product_test():\n' );
  fprintf ( 1, '  Test the value of ||A-U*S*V''||\n' );

  a_norm = norm ( a, 'fro' );

  usv(1:m,1:n) = u(1:m,1:m) * s(1:m,1:n) * v(1:n,1:n)';

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The product U * X * V'':\n' );
  disp ( usv );

  dif_norm = norm ( a - usv, 'fro' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Frobenius Norm of A, A_NORM = %f\n', a_norm );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  ABSOLUTE ERROR for A = U*S*V'':\n' );
  fprintf ( 1, '  Frobenius norm of difference A-U*S*V'' = %f\n', dif_norm );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  RELATIVE ERROR for A = U*S*V'':\n' );
  fprintf ( 1, '  Ratio of DIF_NORM / A_NORM = %f\n', dif_norm / a_norm );

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% timestamp() prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

