GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
svd.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1996-2013 John W. Eaton
4 
5 This file is part of Octave.
6 
7 Octave is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11 
12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
19 <http://www.gnu.org/licenses/>.
20 
21 */
22 
23 #ifdef HAVE_CONFIG_H
24 #include <config.h>
25 #endif
26 
27 #include "CmplxSVD.h"
28 #include "dbleSVD.h"
29 #include "fCmplxSVD.h"
30 #include "floatSVD.h"
31 
32 #include "defun.h"
33 #include "error.h"
34 #include "gripes.h"
35 #include "oct-obj.h"
36 #include "pr-output.h"
37 #include "utils.h"
38 #include "variables.h"
39 
40 static int Vsvd_driver = SVD::GESVD;
41 
42 DEFUN (svd, args, nargout,
43  "-*- texinfo -*-\n\
44 @deftypefn {Built-in Function} {@var{s} =} svd (@var{A})\n\
45 @deftypefnx {Built-in Function} {[@var{U}, @var{S}, @var{V}] =} svd (@var{A})\n\
46 @deftypefnx {Built-in Function} {[@var{U}, @var{S}, @var{V}] =} svd (@var{A}, @var{econ})\n\
47 @cindex singular value decomposition\n\
48 Compute the singular value decomposition of @var{A}\n\
49 @tex\n\
50 $$\n\
51  A = U S V^{\\dagger}\n\
52 $$\n\
53 @end tex\n\
54 @ifnottex\n\
55 \n\
56 @example\n\
57 A = U*S*V'\n\
58 @end example\n\
59 \n\
60 @end ifnottex\n\
61 \n\
62 The function @code{svd} normally returns only the vector of singular values.\n\
63 When called with three return values, it computes\n\
64 @tex\n\
65 $U$, $S$, and $V$.\n\
66 @end tex\n\
67 @ifnottex\n\
68 @var{U}, @var{S}, and @var{V}.\n\
69 @end ifnottex\n\
70 For example,\n\
71 \n\
72 @example\n\
73 svd (hilb (3))\n\
74 @end example\n\
75 \n\
76 @noindent\n\
77 returns\n\
78 \n\
79 @example\n\
80 @group\n\
81 ans =\n\
82 \n\
83  1.4083189\n\
84  0.1223271\n\
85  0.0026873\n\
86 @end group\n\
87 @end example\n\
88 \n\
89 @noindent\n\
90 and\n\
91 \n\
92 @example\n\
93 [u, s, v] = svd (hilb (3))\n\
94 @end example\n\
95 \n\
96 @noindent\n\
97 returns\n\
98 \n\
99 @example\n\
100 @group\n\
101 u =\n\
102 \n\
103  -0.82704 0.54745 0.12766\n\
104  -0.45986 -0.52829 -0.71375\n\
105  -0.32330 -0.64901 0.68867\n\
106 \n\
107 s =\n\
108 \n\
109  1.40832 0.00000 0.00000\n\
110  0.00000 0.12233 0.00000\n\
111  0.00000 0.00000 0.00269\n\
112 \n\
113 v =\n\
114 \n\
115  -0.82704 0.54745 0.12766\n\
116  -0.45986 -0.52829 -0.71375\n\
117  -0.32330 -0.64901 0.68867\n\
118 @end group\n\
119 @end example\n\
120 \n\
121 If given a second argument, @code{svd} returns an economy-sized\n\
122 decomposition, eliminating the unnecessary rows or columns of @var{U} or\n\
123 @var{V}.\n\
124 @seealso{svd_driver, svds, eig, lu, chol, hess, qr, qz}\n\
125 @end deftypefn")
126 {
127  octave_value_list retval;
128 
129  int nargin = args.length ();
130 
131  if (nargin < 1 || nargin > 2 || nargout == 2 || nargout > 3)
132  {
133  print_usage ();
134  return retval;
135  }
136 
137  octave_value arg = args(0);
138 
139  octave_idx_type nr = arg.rows ();
140  octave_idx_type nc = arg.columns ();
141 
142  if (arg.ndims () != 2)
143  {
144  error ("svd: A must be a 2-D matrix");
145  return retval;
146  }
147 
148  bool isfloat = arg.is_single_type ();
149 
150  SVD::type type = ((nargout == 0 || nargout == 1)
152  : (nargin == 2) ? SVD::economy : SVD::std);
153 
154  SVD::driver driver = static_cast<SVD::driver> (Vsvd_driver);
155 
156  if (nr == 0 || nc == 0)
157  {
158  if (isfloat)
159  {
160  switch (type)
161  {
162  case SVD::std:
163  retval(2) = FloatDiagMatrix (nc, nc, 1.0f);
164  retval(1) = FloatMatrix (nr, nc);
165  retval(0) = FloatDiagMatrix (nr, nr, 1.0f);
166  break;
167  case SVD::economy:
168  retval(2) = FloatDiagMatrix (0, nc, 1.0f);
169  retval(1) = FloatMatrix (0, 0);
170  retval(0) = FloatDiagMatrix (nr, 0, 1.0f);
171  break;
172  case SVD::sigma_only: default:
173  retval(0) = FloatMatrix (0, 1);
174  break;
175  }
176  }
177  else
178  {
179  switch (type)
180  {
181  case SVD::std:
182  retval(2) = DiagMatrix (nc, nc, 1.0);
183  retval(1) = Matrix (nr, nc);
184  retval(0) = DiagMatrix (nr, nr, 1.0);
185  break;
186  case SVD::economy:
187  retval(2) = DiagMatrix (0, nc, 1.0);
188  retval(1) = Matrix (0, 0);
189  retval(0) = DiagMatrix (nr, 0, 1.0);
190  break;
191  case SVD::sigma_only: default:
192  retval(0) = Matrix (0, 1);
193  break;
194  }
195  }
196  }
197  else
198  {
199  if (isfloat)
200  {
201  if (arg.is_real_type ())
202  {
203  FloatMatrix tmp = arg.float_matrix_value ();
204 
205  if (! error_state)
206  {
207  if (tmp.any_element_is_inf_or_nan ())
208  {
209  error ("svd: cannot take SVD of matrix containing Inf or NaN values");
210  return retval;
211  }
212 
213  FloatSVD result (tmp, type, driver);
214 
215  FloatDiagMatrix sigma = result.singular_values ();
216 
217  if (nargout == 0 || nargout == 1)
218  {
219  retval(0) = sigma.extract_diag ();
220  }
221  else
222  {
223  retval(2) = result.right_singular_matrix ();
224  retval(1) = sigma;
225  retval(0) = result.left_singular_matrix ();
226  }
227  }
228  }
229  else if (arg.is_complex_type ())
230  {
232 
233  if (! error_state)
234  {
235  if (ctmp.any_element_is_inf_or_nan ())
236  {
237  error ("svd: cannot take SVD of matrix containing Inf or NaN values");
238  return retval;
239  }
240 
241  FloatComplexSVD result (ctmp, type, driver);
242 
243  FloatDiagMatrix sigma = result.singular_values ();
244 
245  if (nargout == 0 || nargout == 1)
246  {
247  retval(0) = sigma.extract_diag ();
248  }
249  else
250  {
251  retval(2) = result.right_singular_matrix ();
252  retval(1) = sigma;
253  retval(0) = result.left_singular_matrix ();
254  }
255  }
256  }
257  }
258  else
259  {
260  if (arg.is_real_type ())
261  {
262  Matrix tmp = arg.matrix_value ();
263 
264  if (! error_state)
265  {
266  if (tmp.any_element_is_inf_or_nan ())
267  {
268  error ("svd: cannot take SVD of matrix containing Inf or NaN values");
269  return retval;
270  }
271 
272  SVD result (tmp, type, driver);
273 
274  DiagMatrix sigma = result.singular_values ();
275 
276  if (nargout == 0 || nargout == 1)
277  {
278  retval(0) = sigma.extract_diag ();
279  }
280  else
281  {
282  retval(2) = result.right_singular_matrix ();
283  retval(1) = sigma;
284  retval(0) = result.left_singular_matrix ();
285  }
286  }
287  }
288  else if (arg.is_complex_type ())
289  {
290  ComplexMatrix ctmp = arg.complex_matrix_value ();
291 
292  if (! error_state)
293  {
294  if (ctmp.any_element_is_inf_or_nan ())
295  {
296  error ("svd: cannot take SVD of matrix containing Inf or NaN values");
297  return retval;
298  }
299 
300  ComplexSVD result (ctmp, type, driver);
301 
302  DiagMatrix sigma = result.singular_values ();
303 
304  if (nargout == 0 || nargout == 1)
305  {
306  retval(0) = sigma.extract_diag ();
307  }
308  else
309  {
310  retval(2) = result.right_singular_matrix ();
311  retval(1) = sigma;
312  retval(0) = result.left_singular_matrix ();
313  }
314  }
315  }
316  else
317  {
318  gripe_wrong_type_arg ("svd", arg);
319  return retval;
320  }
321  }
322  }
323 
324  return retval;
325 }
326 
327 /*
328 %!assert (svd ([1, 2; 2, 1]), [3; 1], sqrt (eps))
329 
330 %!test
331 %! [u, s, v] = svd ([1, 2; 2, 1]);
332 %! x = 1 / sqrt (2);
333 %! assert (u, [-x, -x; -x, x], sqrt (eps));
334 %! assert (s, [3, 0; 0, 1], sqrt (eps));
335 %! assert (v, [-x, x; -x, -x], sqrt (eps));
336 
337 %!test
338 %! a = [1, 2, 3; 4, 5, 6];
339 %! [u, s, v] = svd (a);
340 %! assert (u * s * v', a, sqrt (eps));
341 
342 %!test
343 %! a = [1, 2; 3, 4; 5, 6];
344 %! [u, s, v] = svd (a);
345 %! assert (u * s * v', a, sqrt (eps));
346 
347 %!test
348 %! a = [1, 2, 3; 4, 5, 6];
349 %! [u, s, v] = svd (a, 1);
350 %! assert (u * s * v', a, sqrt (eps));
351 
352 %!test
353 %! a = [1, 2; 3, 4; 5, 6];
354 %! [u, s, v] = svd (a, 1);
355 %! assert (u * s * v', a, sqrt (eps));
356 
357 %!assert (svd (single ([1, 2; 2, 1])), single ([3; 1]), sqrt (eps ("single")))
358 
359 %!test
360 %! [u, s, v] = svd (single ([1, 2; 2, 1]));
361 %! x = single (1 / sqrt (2));
362 %! assert (u, [-x, -x; -x, x], sqrt (eps ("single")));
363 %! assert (s, single ([3, 0; 0, 1]), sqrt (eps ("single")));
364 %! assert (v, [-x, x; -x, -x], sqrt (eps ("single")));
365 
366 %!test
367 %! a = single ([1, 2, 3; 4, 5, 6]);
368 %! [u, s, v] = svd (a);
369 %! assert (u * s * v', a, sqrt (eps ("single")));
370 
371 %!test
372 %! a = single ([1, 2; 3, 4; 5, 6]);
373 %! [u, s, v] = svd (a);
374 %! assert (u * s * v', a, sqrt (eps ("single")));
375 
376 %!test
377 %! a = single ([1, 2, 3; 4, 5, 6]);
378 %! [u, s, v] = svd (a, 1);
379 %! assert (u * s * v', a, sqrt (eps ("single")));
380 
381 %!test
382 %! a = single ([1, 2; 3, 4; 5, 6]);
383 %! [u, s, v] = svd (a, 1);
384 %! assert (u * s * v', a, sqrt (eps ("single")));
385 
386 %!test
387 %! a = zeros (0, 5);
388 %! [u, s, v] = svd (a);
389 %! assert (size (u), [0, 0]);
390 %! assert (size (s), [0, 5]);
391 %! assert (size (v), [5, 5]);
392 
393 %!test
394 %! a = zeros (5, 0);
395 %! [u, s, v] = svd (a, 1);
396 %! assert (size (u), [5, 0]);
397 %! assert (size (s), [0, 0]);
398 %! assert (size (v), [0, 0]);
399 
400 %!error svd ()
401 %!error svd ([1, 2; 4, 5], 2, 3)
402 %!error [u, v] = svd ([1, 2; 3, 4])
403 */
404 
405 DEFUN (svd_driver, args, nargout,
406  "-*- texinfo -*-\n\
407 @deftypefn {Built-in Function} {@var{val} =} svd_driver ()\n\
408 @deftypefnx {Built-in Function} {@var{old_val} =} svd_driver (@var{new_val})\n\
409 @deftypefnx {Built-in Function} {} svd_driver (@var{new_val}, \"local\")\n\
410 Query or set the underlying @sc{lapack} driver used by @code{svd}.\n\
411 Currently recognized values are @qcode{\"gesvd\"} and @qcode{\"gesdd\"}. \n\
412 The default is @qcode{\"gesvd\"}.\n\
413 \n\
414 When called from inside a function with the @qcode{\"local\"} option, the\n\
415 variable is changed locally for the function and any subroutines it calls. \n\
416 The original variable value is restored when exiting the function.\n\
417 @seealso{svd}\n\
418 @end deftypefn")
419 {
420  static const char *driver_names[] = { "gesvd", "gesdd", 0 };
421 
422  return SET_INTERNAL_VARIABLE_CHOICES (svd_driver, driver_names);
423 }
424 
425 /*
426 %!test
427 %! A = [1+1i, 1-1i, 0; 0, 2, 0; 1i, 1i, 1+2i];
428 %! old_driver = svd_driver ("gesvd");
429 %! [U1, S1, V1] = svd (A);
430 %! svd_driver ("gesdd");
431 %! [U2, S2, V2] = svd (A);
432 %! assert (U1, U2, 5*eps);
433 %! assert (S1, S2, 5*eps);
434 %! assert (V1, V2, 5*eps);
435 %! svd_driver (old_driver);
436 */