Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

svd: compute numerically when matrix has Floats and Integers #766

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

svd: compute numerically when matrix has Floats and Integers #766

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
185 changes: 174 additions & 11 deletions inst/@sym/svd.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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.
%%
Expand All @@ -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]);
Expand All @@ -37,35 +51,159 @@
%%
%% @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)'
'if svecs:'
' (U, S, Vt) = tmp'
' U = Matrix(U.rows, U.cols, lambda i,j: U[i, j])'
' 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];
Expand Down Expand Up @@ -99,3 +237,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])