octave  3.8.2
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-3.8.2.tar.gz  ("inofficial" and yet experimental doxygen-generated source code documentation)  

 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-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 "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}. Singular values less than\n\
44 @var{tol} are ignored.\n\
45 \n\
46 If the second argument is omitted, it is taken to be\n\
47 \n\
48 @example\n\
49 tol = max (size (@var{x})) * sigma_max (@var{x}) * eps,\n\
50 @end example\n\
51 \n\
52 @noindent\n\
53 where @code{sigma_max (@var{x})} is the maximal singular value of @var{x}.\n\
54 @end deftypefn")
55 {
57 
58  int nargin = args.length ();
59 
60  if (nargin < 1 || nargin > 2)
61  {
62  print_usage ();
63  return retval;
64  }
65 
66  octave_value arg = args(0);
67 
68  int arg_is_empty = empty_arg ("pinv", arg.rows (), arg.columns ());
69 
70  if (arg_is_empty < 0)
71  return retval;
72  else if (arg_is_empty > 0)
73  return octave_value (Matrix ());
74 
75  bool isfloat = arg.is_single_type ();
76 
77  if (arg.is_diag_matrix ())
78  {
79  if (nargin == 2)
80  warning ("pinv: tol is ignored for diagonal matrices");
81 
82  if (arg.is_complex_type ())
83  {
84  if (isfloat)
86  else
87  retval = arg.complex_diag_matrix_value ().pseudo_inverse ();
88  }
89  else
90  {
91  if (isfloat)
92  retval = arg.float_diag_matrix_value ().pseudo_inverse ();
93  else
94  retval = arg.diag_matrix_value ().pseudo_inverse ();
95  }
96  }
97  else if (arg.is_perm_matrix ())
98  {
99  retval = arg.perm_matrix_value ().inverse ();
100  }
101  else if (isfloat)
102  {
103  float tol = 0.0;
104  if (nargin == 2)
105  tol = args(1).float_value ();
106 
107  if (error_state)
108  return retval;
109 
110  if (tol < 0.0)
111  {
112  error ("pinv: TOL must be greater than zero");
113  return retval;
114  }
115 
116  if (arg.is_real_type ())
117  {
119 
120  if (! error_state)
121  retval = m.pseudo_inverse (tol);
122  }
123  else if (arg.is_complex_type ())
124  {
126 
127  if (! error_state)
128  retval = m.pseudo_inverse (tol);
129  }
130  else
131  {
132  gripe_wrong_type_arg ("pinv", arg);
133  }
134  }
135  else
136  {
137  double tol = 0.0;
138  if (nargin == 2)
139  tol = args(1).double_value ();
140 
141  if (error_state)
142  return retval;
143 
144  if (tol < 0.0)
145  {
146  error ("pinv: TOL must be greater than zero");
147  return retval;
148  }
149 
150  if (arg.is_real_type ())
151  {
152  Matrix m = arg.matrix_value ();
153 
154  if (! error_state)
155  retval = m.pseudo_inverse (tol);
156  }
157  else if (arg.is_complex_type ())
158  {
160 
161  if (! error_state)
162  retval = m.pseudo_inverse (tol);
163  }
164  else
165  {
166  gripe_wrong_type_arg ("pinv", arg);
167  }
168  }
169 
170  return retval;
171 }
172 
173 /*
174 %!shared a, b, tol, hitol, d, u, x, y
175 %! a = reshape (rand*[1:16], 4, 4); # Rank 2 matrix
176 %! b = pinv (a);
177 %! tol = 4e-14;
178 %! hitol = 40*sqrt (eps);
179 %! d = diag ([rand, rand, hitol, hitol]);
180 %! u = rand (4); # Could be singular by freak accident
181 %! x = inv (u)*d*u;
182 %! y = pinv (x, sqrt (eps));
183 %!
184 %!assert (a*b*a, a, tol)
185 %!assert (b*a*b, b, tol)
186 %!assert ((b*a)', b*a, tol)
187 %!assert ((a*b)', a*b, tol)
188 %!assert (x*y*x, x, -hitol)
189 %!assert (y*x*y, y, -hitol)
190 %!assert ((x*y)', x*y, hitol)
191 %!assert ((y*x)', y*x, hitol)
192 */
DiagMatrix pseudo_inverse(void) const
Definition: dDiagMatrix.cc:295
FloatComplexDiagMatrix float_complex_diag_matrix_value(bool force=false) const
Definition: ov.h:830
bool is_real_type(void) const
Definition: ov.h:641
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:466
FloatComplexMatrix float_complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:785
endif else print_usage
Definition: nargoutchk.m:102
args
Definition: mk_doc_cache.m:17
endif single tol
Definition: cplxpair.m:66
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}. Singular values less than\n\ @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
else arg
Definition: clabel.m:75
bool is_perm_matrix(void) const
Definition: ov.h:552
function retval
Definition: mput.m:29
int empty_arg(const char *, octave_idx_type nr, octave_idx_type nc)
Definition: utils.cc:245
ComplexDiagMatrix complex_diag_matrix_value(bool force=false) const
Definition: ov.h:826
FloatComplexMatrix pseudo_inverse(float tol=0.0) const
Definition: fCMatrix.cc:1194
octave_idx_type columns(void) const
Definition: ov.h:468
FloatComplexDiagMatrix pseudo_inverse(void) const
PermMatrix inverse(void) const
Definition: PermMatrix.cc:95
FloatDiagMatrix float_diag_matrix_value(bool force=false) const
Definition: ov.h:823
int error_state
Definition: error.cc:100
bool is_complex_type(void) const
Definition: ov.h:644
unwind_protect warning("off")
Definition: fail.m:145
Definition: dMatrix.h:34
Matrix matrix_value(bool frc_str_conv=false) const
Definition: ov.h:763
return octave_value()
endif nargin
Definition: daspect.m:63
DiagMatrix diag_matrix_value(bool force=false) const
Definition: ov.h:820
ComplexMatrix pseudo_inverse(double tol=0.0) const
Definition: CMatrix.cc:1186
ComplexDiagMatrix pseudo_inverse(void) const
Definition: CDiagMatrix.cc:386
ComplexMatrix complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:781
FloatMatrix float_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:766
FloatDiagMatrix pseudo_inverse(void) const
Definition: fDiagMatrix.cc:295
PermMatrix perm_matrix_value(void) const
Definition: ov.h:833
endif if(nargin==1||(nargin > 1 &&islogical(cond)&&ischar(varargin{1}))) if((!isnumeric(cond)&&!islogical(cond))||!all(cond( error("assert %s failed", argin)
Definition: assert.m:79
bool is_single_type(void) const
Definition: ov.h:601
bool is_diag_matrix(void) const
Definition: ov.h:549
Matrix pseudo_inverse(double tol=0.0) const
Definition: dMatrix.cc:854
FloatMatrix pseudo_inverse(float tol=0.0) const
Definition: fMatrix.cc:858