svd: compute numerically when matrix has Floats and Integers

Maybe we want to do this upstream in SymPy, so sort of RFC for now.

inst/@sym/svd.m

@@ -1,4 +1,4 @@
-%% Copyright (C) 2014, 2016 Colin B. Macdonald
+%% Copyright (C) 2014, 2016-2017 Colin B. Macdonald
 %%
 %% This file is part of OctSymPy.
 %%
@@ -20,11 +20,25 @@
 %% @documentencoding UTF-8
 %% @deftypemethod  @@sym {@var{S} =} svd (@var{A})
 %% @deftypemethodx @@sym {[@var{U}, @var{S}, @var{V}] =} svd (@var{A})
+%% @deftypemethodx @@sym {[@var{U}, @var{S}, @var{V}] =} svd (@var{A}, 'econ')
 %% Symbolic singular value decomposition.
 %%
-%% The SVD: U*S*V' = A
+%% The SVD is the factorization of matrix
+%% @iftex
+%% @math{A} into @math{U S V^T = A},
+%% where @math{S} is a diagonal matrix of @emph{singular values},
+%% and @math{U} and @math{V}
+%% @end iftex
+%% @ifnottex
+%% A into U*S*V' = A,
+%% where S is a diagonal matrix of @emph{singular values},
+%% and U and V
+%% @end ifnottex
+%% are orthogonal matrices whose columns form the left and right
+%% @emph{singular vectors}.
 %%
-%% Singular values example:
+%% When the matrix contains symbols, expressions, rational numbers, and other
+%% things, this command finds the singular values via symbolic manipulation:
 %% @example
 %% @group
 %% A = sym([1 0; 3 0]);
@@ -37,35 +51,161 @@
 %%
 %% @end group
 %% @end example
+%% Currently only the singular values (not the singular vectors) are supported
+%% in symbolic mode.
 %%
-%% FIXME: currently only singular values, not singular vectors.
-%% Should add full SVD to sympy.
 %%
+%% If the matrix contains Float entries (@pxref{vpa}) (and possibly Integers),
+%% the SVD is computing numerically in variable precision
+%% arithmetic in the precision given by @pxref{digits}.
+%% The singular values and singular vectors can be computed in this mode.
+%% Example:
+%% @example
+%% @group
+%% A = vpa (3*hilb (sym(3)));
+%% [U, S, V] = svd (A)
+%%   @result{} U = (sym 3×3 matrix)
+%%       ...
+%%   @result{} S = (sym 3×3 matrix)
+%%       ...
+%%   @result{} V = (sym 3×3 matrix)
+%%       ...
+%% @end group
+%%
+%% @group
+%% diag(S)
+%%   @result{} (sym 3×1 matrix)
+%%       ⎡ 4.2249567813709618726124675083017 ⎤
+%%       ⎢                                   ⎥
+%%       ⎢ 0.3669811975617175396944034278698 ⎥
+%%       ⎢                                   ⎥
+%%       ⎣0.008062021067320587693129063828546⎦
+%% @end group
+%% @end example
+%%
+%% Next, extract one singular value and associated left/right
+%% singular vectors:
+%% @example
+%% @group
+%% sv = S(1, 1)
+%% u = U(:, 1)
+%% v = V(:, 1)
+%%   @result{} sv = (sym) 4.2249567813709618726124675083017
+%%   @result{} u = (sym 3×1 matrix)
+%%       ⎡-0.82704492697200940922027703647284⎤
+%%       ⎢                                   ⎥
+%%       ⎢-0.45986390436554392104852568981886⎥
+%%       ⎢                                   ⎥
+%%       ⎣-0.32329843524449897629157179151973⎦
+%%   @result{} v = (sym 3×1 matrix)
+%%       ⎡-0.82704492697200940922027703647284⎤
+%%       ⎢                                   ⎥
+%%       ⎢-0.45986390436554392104852568981886⎥
+%%       ⎢                                   ⎥
+%%       ⎣-0.32329843524449897629157179151973⎦
+%% @end group
+%% @end example
+%%
+%% Check the SVD is satisfied to high-precision:
+%% @example
+%% @group
+%% sv*u - A*v
+%%   @result{} (sym 3×1 matrix)
+%%       ⎡-9.2444637330587320946686941244077e-33⎤
+%%       ⎢                                      ⎥
+%%       ⎢-3.0814879110195773648895647081359e-33⎥
+%%       ⎢                                      ⎥
+%%       ⎣-3.0814879110195773648895647081359e-33⎦
+%% @end group
+%% @end example
+%%
+%% If the @qcode{'econ'} keyboard is passed, an ``economy size''
+%% SVD is returned (@pxref{svd}).
 %% @seealso{svd, @@sym/eig}
 %% @end deftypemethod
 
 
-function [S, varargout] = svd(A)
+function [S, varargout] = svd(A, econ)
 
-  if (nargin >= 2)
-    error('svd: economy-size not supported yet')
+  if (nargin == 1)
+    econ = false;
+  elseif (nargin == 2)
+    if (isnumeric(econ) && econ == 0)
+      error('svd: auto econ mode ("0") is not yet supported')
+    else
+      econ = true;
+    end
+  else
+    print_usage ();
   end
 
-  if (nargout >= 2)
-    error('svd: singular vectors not yet computed by sympy')
+  if (nargout <= 1)
+    svecs = false;
+  elseif (nargout == 3)
+    svecs = true;
+  else
+    print_usage ();
   end
 
+
+  cmd = { 'A, = _ins'
+          'A = A if A.is_Matrix else Matrix([A])'
+          'return (any([x.is_Float for x in A]) and'
+          '       all([x.is_Float or x.is_Integer for x in A]))' };
+  is_vpa_matrix = python_cmd (cmd, sym(A));
+
+  if (is_vpa_matrix)
+    myd = digits ();   % TODO: or take from the object itself
+    cmd = { '(A, svecs, econ, digits) = _ins'
+            'A = A if A.is_Matrix else Matrix([A])'
+            'import mpmath'
+            'mpmath.mp.dps = digits'
+            'tmp = mpmath.svd(mpmath.matrix(A), full_matrices=(not econ), compute_uv=svecs)'
+            '#dbout(tmp)'
+            'if svecs:'
+            '    (U, S, Vt) = tmp'
+            '    U = Matrix(U.rows, U.cols, lambda i,j: U[i, j])'
+            '    #assert U.is_Matrix'  % made from a copy of A
+            '    S = Matrix(S)'
+            '    m, n = A.shape'
+            '    r, c = (m, n) if not econ else (min(m,n),)*2'
+            '    S = Matrix(r, c, lambda i,j: S[i] if i == j else 0)'
+            '    V = Vt.transpose()'  % TODO: or transpose_conj?
+            '    V = Matrix(V.rows, V.cols, lambda i,j: V[i, j])'
+            'else:'
+            '    S = Matrix(tmp)'
+            '    U, V = None, None'
+            'return (U, S, V)' };
+    [U, S, V] = python_cmd (cmd, sym(A), svecs, econ, myd);
+
+    if (nargout >= 2)
+      varargout{1} = S;
+      varargout{2} = V;
+      S = U;
+    end
+
+  else
+    if (svecs)
+      error ('svd: singular vectors not yet implemented for non-vpa matrices')
+    end
+    if (econ)
+      error ('svd: "economy size" not yet implemented for non-vpa matrices')
+    end
+
   cmd = { '(A,) = _ins'
           'if not A.is_Matrix:'
           '    A = sp.Matrix([A])'
           'L = sp.Matrix(A.singular_values())'
           'return L,' };
 
   S = python_cmd (cmd, sym(A));
-
+  end
 end
 
 
+%!error <Invalid> svd (sym(1), 2, 3)
+%!error <Invalid> [a, b] = svd (sym(1))
+
 %!test
 %! % basic
 %! A = [1 2; 3 4];
@@ -99,3 +239,28 @@
 %%! A = [x 0; sym(0) 2*x]
 %%! [u,s,v] = cond(A)
 %%! assert (false)
+
+%!test
+%! % econ & non-square matrices
+%! A = vpa([1 2 4; 1 2 4]);
+%! S = svd (A);
+%! assert (size (S), [2 1])
+%! [U, S, V] = svd (A)
+%! assert (size (U), [2 2])
+%! assert (size (S), [2 3])
+%! assert (size (V), [3 3])
+%! [U, S, V] = svd (A, 'econ');
+%! assert (size (U), [2 2])
+%! assert (size (S), [2 2])
+%! assert (size (V), [3 2])
+%! A = A';
+%! S = svd (A);
+%! assert (size (S), [2 1])
+%! [U, S, V] = svd (A, 'econ');
+%! assert (size (U), [3 2])
+%! assert (size (S), [2 2])
+%! assert (size (V), [2 2])
+%! [U, S, V] = svd (A);
+%! assert (size (U), [3 3])
+%! assert (size (S), [3 2])
+%! assert (size (V), [2 2])