octave  4.0.0
About: Octave is a high-level language for linear and nonlinear numerical computations mostly compatible with Matlab (note: Octave 3.x is significantly different from Octave 2.1).
  Fossies Dox: octave-4.0.0.tar.gz  ("inofficial" and yet experimental doxygen-generated source code documentation)  

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 {
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  {
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
endif else print_usage
Definition: nargoutchk.m:103
DEFUN(pinv, args,,"-*- texinfo -*-\n\ @deftypefn {Built-in Function} {} pinv (@var{x})\n\ @deftypefnx {Built-in Function} {} pinv (@var{x}, @var{tol})\n\ Return the pseudoinverse of @var{x}.\n\ \n\ Singular values less than @var{tol} are ignored.\n\ \n\ If the second argument is omitted, it is taken to be\n\ \n\ @example\n\ tol = max (size (@var{x})) * sigma_max (@var{x}) * eps,\n\ @end example\n\ \n\ @noindent\n\ where @code{sigma_max (@var{x})} is the maximal singular value of @var{x}.\n\ @end deftypefn")
Definition: pinv.cc:39
nargin
Definition: annotation.m:152
endif single tol
Definition: cplxpair.m:68
else arg
Definition: clabel.m:75
bool is_perm_matrix(void) const
Definition: ov.h:559
function retval
Definition: mput.m:30
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
args
Definition: genpropdoc.m:1174
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
Definition: dMatrix.h:35
Matrix matrix_value(bool frc_str_conv=false) const
Definition: ov.h:773
return octave_value()
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
endif if(nargin==1||(nargin > 1 &&islogical(cond)&&ischar(varargin{1}))) if((!isnumeric(cond)&&!islogical(cond))||isempty(cond)||!all(cond( error("assert %s failed", argin)
Definition: assert.m:92
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
FloatMatrix pseudo_inverse(float tol=0.0) const
Definition: fMatrix.cc:877