GNU Octave  4.0.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
pinv.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1996-2015 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 "defun.h"
28 #include "error.h"
29 #include "gripes.h"
30 #include "oct-obj.h"
31 #include "utils.h"
32 #include "ops.h"
33 #include "ov-re-diag.h"
34 #include "ov-cx-diag.h"
35 #include "ov-flt-re-diag.h"
36 #include "ov-flt-cx-diag.h"
37 #include "ov-perm.h"
38 
39 DEFUN (pinv, args, ,
40  "-*- texinfo -*-\n\
41 @deftypefn {Built-in Function} {} pinv (@var{x})\n\
42 @deftypefnx {Built-in Function} {} pinv (@var{x}, @var{tol})\n\
43 Return the pseudoinverse of @var{x}.\n\
44 \n\
45 Singular values less than @var{tol} are ignored.\n\
46 \n\
47 If the second argument is omitted, it is taken to be\n\
48 \n\
49 @example\n\
50 tol = max (size (@var{x})) * sigma_max (@var{x}) * eps,\n\
51 @end example\n\
52 \n\
53 @noindent\n\
54 where @code{sigma_max (@var{x})} is the maximal singular value of @var{x}.\n\
55 @end deftypefn")
56 {
57  octave_value retval;
58 
59  int nargin = args.length ();
60 
61  if (nargin < 1 || nargin > 2)
62  {
63  print_usage ();
64  return retval;
65  }
66 
67  octave_value arg = args(0);
68 
69  int arg_is_empty = empty_arg ("pinv", arg.rows (), arg.columns ());
70 
71  if (arg_is_empty < 0)
72  return retval;
73  else if (arg_is_empty > 0)
74  return octave_value (Matrix ());
75 
76  bool isfloat = arg.is_single_type ();
77 
78  if (arg.is_diag_matrix ())
79  {
80  if (isfloat)
81  {
82  float tol = 0.0;
83  if (nargin == 2)
84  tol = args(1).float_value ();
85 
86  if (error_state)
87  return retval;
88 
89  if (tol < 0.0)
90  {
91  error ("pinv: TOL must be greater than zero");
92  return retval;
93  }
94 
95  if (arg.is_real_type ())
96  retval = arg.float_diag_matrix_value ().pseudo_inverse (tol);
97  else
98  retval = arg.float_complex_diag_matrix_value ().pseudo_inverse (tol);
99  }
100  else
101  {
102  double tol = 0.0;
103  if (nargin == 2)
104  tol = args(1).double_value ();
105 
106  if (error_state)
107  return retval;
108 
109  if (tol < 0.0)
110  {
111  error ("pinv: TOL must be greater than zero");
112  return retval;
113  }
114 
115  if (arg.is_real_type ())
116  retval = arg.diag_matrix_value ().pseudo_inverse (tol);
117  else
118  retval = arg.complex_diag_matrix_value ().pseudo_inverse (tol);
119  }
120  }
121  else if (arg.is_perm_matrix ())
122  {
123  retval = arg.perm_matrix_value ().inverse ();
124  }
125  else if (isfloat)
126  {
127  float tol = 0.0;
128  if (nargin == 2)
129  tol = args(1).float_value ();
130 
131  if (error_state)
132  return retval;
133 
134  if (tol < 0.0)
135  {
136  error ("pinv: TOL must be greater than zero");
137  return retval;
138  }
139 
140  if (arg.is_real_type ())
141  {
142  FloatMatrix m = arg.float_matrix_value ();
143 
144  if (! error_state)
145  retval = m.pseudo_inverse (tol);
146  }
147  else if (arg.is_complex_type ())
148  {
150 
151  if (! error_state)
152  retval = m.pseudo_inverse (tol);
153  }
154  else
155  {
156  gripe_wrong_type_arg ("pinv", arg);
157  }
158  }
159  else
160  {
161  double tol = 0.0;
162  if (nargin == 2)
163  tol = args(1).double_value ();
164 
165  if (error_state)
166  return retval;
167 
168  if (tol < 0.0)
169  {
170  error ("pinv: TOL must be greater than zero");
171  return retval;
172  }
173 
174  if (arg.is_real_type ())
175  {
176  Matrix m = arg.matrix_value ();
177 
178  if (! error_state)
179  retval = m.pseudo_inverse (tol);
180  }
181  else if (arg.is_complex_type ())
182  {
184 
185  if (! error_state)
186  retval = m.pseudo_inverse (tol);
187  }
188  else
189  {
190  gripe_wrong_type_arg ("pinv", arg);
191  }
192  }
193 
194  return retval;
195 }
196 
197 /*
198 %!shared a, b, tol, hitol, d, u, x, y
199 %! a = reshape (rand*[1:16], 4, 4); # Rank 2 matrix
200 %! b = pinv (a);
201 %! tol = 4e-14;
202 %! hitol = 40*sqrt (eps);
203 %! d = diag ([rand, rand, hitol, hitol]);
204 %! u = rand (4); # Could be singular by freak accident
205 %! x = inv (u)*d*u;
206 %! y = pinv (x, sqrt (eps));
207 %!
208 %!assert (a*b*a, a, tol)
209 %!assert (b*a*b, b, tol)
210 %!assert ((b*a)', b*a, tol)
211 %!assert ((a*b)', a*b, tol)
212 %!assert (x*y*x, x, -hitol)
213 %!assert (y*x*y, y, -hitol)
214 %!assert ((x*y)', x*y, hitol)
215 %!assert ((y*x)', y*x, hitol)
216 
217 ## Clear shared variables
218 %!shared
219 
220 ## Test pinv for Diagonal matrices
221 %!test
222 %! x = diag ([3 2 1 0 -0.5]);
223 %! y = pinv (x);
224 %! assert (typeinfo (y)(1:8), "diagonal");
225 %! assert (isa (y, "double"));
226 %! assert (diag (y), [1/3, 1/2, 1, 0 1/-0.5]');
227 %! y = pinv (x, 1);
228 %! assert (diag (y), [1/3 1/2 1 0 0]');
229 %! y = pinv (x, 2);
230 %! assert (diag (y), [1/3 1/2 0 0 0]');
231 
232 */
FloatComplexDiagMatrix float_complex_diag_matrix_value(bool force=false) const
Definition: ov.h:840
bool is_real_type(void) const
Definition: ov.h:651
void gripe_wrong_type_arg(const char *name, const char *s, bool is_error)
Definition: gripes.cc:135
octave_idx_type rows(void) const
Definition: ov.h:473
FloatComplexMatrix float_complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:795
OCTINTERP_API void print_usage(void)
Definition: defun.cc:51
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:44
void error(const char *fmt,...)
Definition: error.cc:476
bool is_perm_matrix(void) const
Definition: ov.h:559
DiagMatrix pseudo_inverse(double tol=0.0) const
Definition: dDiagMatrix.cc:295
int empty_arg(const char *, octave_idx_type nr, octave_idx_type nc)
Definition: utils.cc:246
ComplexDiagMatrix pseudo_inverse(double tol=0.0) const
Definition: CDiagMatrix.cc:386
ComplexDiagMatrix complex_diag_matrix_value(bool force=false) const
Definition: ov.h:836
FloatComplexMatrix pseudo_inverse(float tol=0.0) const
Definition: fCMatrix.cc:1227
octave_idx_type columns(void) const
Definition: ov.h:475
FloatComplexDiagMatrix pseudo_inverse(float tol=0.0f) const
PermMatrix inverse(void) const
Definition: PermMatrix.cc:132
FloatDiagMatrix float_diag_matrix_value(bool force=false) const
Definition: ov.h:833
int error_state
Definition: error.cc:101
bool is_complex_type(void) const
Definition: ov.h:654
octave_idx_type length(void) const
Definition: ov.cc:1525
Definition: dMatrix.h:35
Matrix matrix_value(bool frc_str_conv=false) const
Definition: ov.h:773
double arg(double x)
Definition: lo-mappers.h:37
DiagMatrix diag_matrix_value(bool force=false) const
Definition: ov.h:830
FloatDiagMatrix pseudo_inverse(float tol=0.0f) const
Definition: fDiagMatrix.cc:295
ComplexMatrix pseudo_inverse(double tol=0.0) const
Definition: CMatrix.cc:1220
ComplexMatrix complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:791
FloatMatrix float_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:776
PermMatrix perm_matrix_value(void) const
Definition: ov.h:843
bool is_single_type(void) const
Definition: ov.h:611
bool is_diag_matrix(void) const
Definition: ov.h:556
Matrix pseudo_inverse(double tol=0.0) const
Definition: dMatrix.cc:870
return octave_value(v1.char_array_value().concat(v2.char_array_value(), ra_idx),((a1.is_sq_string()||a2.is_sq_string())? '\'': '"'))
FloatMatrix pseudo_inverse(float tol=0.0) const
Definition: fMatrix.cc:877