Skip to content

Commit 7061471

Browse files
cbm755cbm755
committed
svd: compute numerically when matrix has Floats and Integers
Maybe we want to do this upstream in SymPy, so sort of RFC for now.
1 parent fca5944 commit 7061471

File tree

1 file changed

+174
-11
lines changed

1 file changed

+174
-11
lines changed

‎inst/@sym/svd.m‎

Lines changed: 174 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
%% Copyright (C) 2014, 2016 Colin B. Macdonald
1+
%% Copyright (C) 2014, 2016-2017 Colin B. Macdonald
22
%%
33
%% This file is part of OctSymPy.
44
%%
@@ -20,11 +20,25 @@
2020
%% @documentencoding UTF-8
2121
%% @deftypemethod @@sym {@var{S} =} svd (@var{A})
2222
%% @deftypemethodx @@sym {[@var{U}, @var{S}, @var{V}] =} svd (@var{A})
23+
%% @deftypemethodx @@sym {[@var{U}, @var{S}, @var{V}] =} svd (@var{A}, 'econ')
2324
%% Symbolic singular value decomposition.
2425
%%
25-
%% The SVD: U*S*V' = A
26+
%% The SVD is the factorization of matrix
27+
%% @iftex
28+
%% @math{A} into @math{U S V^T = A},
29+
%% where @math{S} is a diagonal matrix of @emph{singular values},
30+
%% and @math{U} and @math{V}
31+
%% @end iftex
32+
%% @ifnottex
33+
%% A into U*S*V' = A,
34+
%% where S is a diagonal matrix of @emph{singular values},
35+
%% and U and V
36+
%% @end ifnottex
37+
%% are orthogonal matrices whose columns form the left and right
38+
%% @emph{singular vectors}.
2639
%%
27-
%% Singular values example:
40+
%% When the matrix contains symbols, expressions, rational numbers, and other
41+
%% things, this command finds the singular values via symbolic manipulation:
2842
%% @example
2943
%% @group
3044
%% A = sym([1 0; 3 0]);
@@ -37,35 +51,159 @@
3751
%%
3852
%% @end group
3953
%% @end example
54+
%% Currently only the singular values (not the singular vectors) are supported
55+
%% in symbolic mode.
4056
%%
41-
%% FIXME: currently only singular values, not singular vectors.
42-
%% Should add full SVD to sympy.
4357
%%
58+
%% If the matrix contains Float entries (@pxref{vpa}) (and possibly Integers),
59+
%% the SVD is computing numerically in variable precision
60+
%% arithmetic in the precision given by @pxref{digits}.
61+
%% The singular values and singular vectors can be computed in this mode.
62+
%% Example:
63+
%% @example
64+
%% @group
65+
%% A = vpa (3*hilb (sym(3)));
66+
%% [U, S, V] = svd (A)
67+
%% @result{} U = (sym 3×3 matrix)
68+
%% ...
69+
%% @result{} S = (sym 3×3 matrix)
70+
%% ...
71+
%% @result{} V = (sym 3×3 matrix)
72+
%% ...
73+
%% @end group
74+
%%
75+
%% @group
76+
%% diag(S)
77+
%% @result{} (sym 3×1 matrix)
78+
%% ⎡ 4.2249567813709618726124675083017 ⎤
79+
%% ⎢ ⎥
80+
%% ⎢ 0.3669811975617175396944034278698 ⎥
81+
%% ⎢ ⎥
82+
%% ⎣0.008062021067320587693129063828546⎦
83+
%% @end group
84+
%% @end example
85+
%%
86+
%% Next, extract one singular value and associated left/right
87+
%% singular vectors:
88+
%% @example
89+
%% @group
90+
%% sv = S(1, 1)
91+
%% u = U(:, 1)
92+
%% v = V(:, 1)
93+
%% @result{} sv = (sym) 4.2249567813709618726124675083017
94+
%% @result{} u = (sym 3×1 matrix)
95+
%% ⎡-0.82704492697200940922027703647284⎤
96+
%% ⎢ ⎥
97+
%% ⎢-0.45986390436554392104852568981886⎥
98+
%% ⎢ ⎥
99+
%% ⎣-0.32329843524449897629157179151973⎦
100+
%% @result{} v = (sym 3×1 matrix)
101+
%% ⎡-0.82704492697200940922027703647284⎤
102+
%% ⎢ ⎥
103+
%% ⎢-0.45986390436554392104852568981886⎥
104+
%% ⎢ ⎥
105+
%% ⎣-0.32329843524449897629157179151973⎦
106+
%% @end group
107+
%% @end example
108+
%%
109+
%% Check the SVD is satisfied to high-precision:
110+
%% @example
111+
%% @group
112+
%% sv*u - A*v
113+
%% @result{} (sym 3×1 matrix)
114+
%% ⎡-9.2444637330587320946686941244077e-33⎤
115+
%% ⎢ ⎥
116+
%% ⎢-3.0814879110195773648895647081359e-33⎥
117+
%% ⎢ ⎥
118+
%% ⎣-3.0814879110195773648895647081359e-33⎦
119+
%% @end group
120+
%% @end example
121+
%%
122+
%% If the @qcode{'econ'} keyboard is passed, an ``economy size''
123+
%% SVD is returned (@pxref{svd}).
44124
%% @seealso{svd, @@sym/eig}
45125
%% @end deftypemethod
46126

47127

48-
function [S, varargout] = svd(A)
128+
function [S, varargout] = svd(A, econ)
49129

50-
if (nargin >= 2)
51-
error('svd: economy-size not supported yet')
130+
if (nargin == 1)
131+
econ = false;
132+
elseif (nargin == 2)
133+
if (isnumeric(econ) && econ == 0)
134+
error('svd: auto econ mode ("0") is not yet supported')
135+
else
136+
econ = true;
137+
end
138+
else
139+
print_usage ();
52140
end
53141

54-
if (nargout >= 2)
55-
error('svd: singular vectors not yet computed by sympy')
142+
if (nargout <= 1)
143+
svecs = false;
144+
elseif (nargout == 3)
145+
svecs = true;
146+
else
147+
print_usage ();
56148
end
57149

150+
151+
cmd = { 'A, = _ins'
152+
'A = A if A.is_Matrix else Matrix([A])'
153+
'return (any([x.is_Float for x in A]) and'
154+
' all([x.is_Float or x.is_Integer for x in A]))' };
155+
is_vpa_matrix = python_cmd (cmd, sym(A));
156+
157+
if (is_vpa_matrix)
158+
myd = digits (); % TODO: or take from the object itself
159+
cmd = { '(A, svecs, econ, digits) = _ins'
160+
'A = A if A.is_Matrix else Matrix([A])'
161+
'import mpmath'
162+
'mpmath.mp.dps = digits'
163+
'tmp = mpmath.svd(mpmath.matrix(A), full_matrices=(not econ), compute_uv=svecs)'
164+
'if svecs:'
165+
' (U, S, Vt) = tmp'
166+
' U = Matrix(U.rows, U.cols, lambda i,j: U[i, j])'
167+
' S = Matrix(S)'
168+
' m, n = A.shape'
169+
' r, c = (m, n) if not econ else (min(m,n),)*2'
170+
' S = Matrix(r, c, lambda i,j: S[i] if i == j else 0)'
171+
' V = Vt.transpose()' % TODO: or transpose_conj?
172+
' V = Matrix(V.rows, V.cols, lambda i,j: V[i, j])'
173+
'else:'
174+
' S = Matrix(tmp)'
175+
' U, V = None, None'
176+
'return (U, S, V)' };
177+
[U, S, V] = python_cmd (cmd, sym(A), svecs, econ, myd);
178+
179+
if (nargout >= 2)
180+
varargout{1} = S;
181+
varargout{2} = V;
182+
S = U;
183+
end
184+
185+
else
186+
if (svecs)
187+
error ('svd: singular vectors not yet implemented for non-vpa matrices')
188+
end
189+
if (econ)
190+
error ('svd: "economy size" not yet implemented for non-vpa matrices')
191+
end
192+
58193
cmd = { '(A,) = _ins'
59194
'if not A.is_Matrix:'
60195
' A = sp.Matrix([A])'
61196
'L = sp.Matrix(A.singular_values())'
62197
'return L,' };
63198

64199
S = python_cmd (cmd, sym(A));
65-
200+
end
66201
end
67202

68203

204+
%!error <Invalid> svd (sym(1), 2, 3)
205+
%!error <Invalid> [a, b] = svd (sym(1))
206+
69207
%!test
70208
%! % basic
71209
%! A = [1 2; 3 4];
@@ -99,3 +237,28 @@
99237
%%! A = [x 0; sym(0) 2*x]
100238
%%! [u,s,v] = cond(A)
101239
%%! assert (false)
240+
241+
%!test
242+
%! % econ & non-square matrices
243+
%! A = vpa([1 2 4; 1 2 4]);
244+
%! S = svd (A);
245+
%! assert (size (S), [2 1])
246+
%! [U, S, V] = svd (A);
247+
%! assert (size (U), [2 2])
248+
%! assert (size (S), [2 3])
249+
%! assert (size (V), [3 3])
250+
%! [U, S, V] = svd (A, 'econ');
251+
%! assert (size (U), [2 2])
252+
%! assert (size (S), [2 2])
253+
%! assert (size (V), [3 2])
254+
%! A = A';
255+
%! S = svd (A);
256+
%! assert (size (S), [2 1])
257+
%! [U, S, V] = svd (A, 'econ');
258+
%! assert (size (U), [3 2])
259+
%! assert (size (S), [2 2])
260+
%! assert (size (V), [2 2])
261+
%! [U, S, V] = svd (A);
262+
%! assert (size (U), [3 3])
263+
%! assert (size (S), [3 2])
264+
%! assert (size (V), [2 2])

0 commit comments

Comments
 (0)