1 //////////////////////////////////////////////////////////////////////// 2 // 3 // Copyright (C) 1996-2021 The Octave Project Developers 4 // 5 // See the file COPYRIGHT.md in the top-level directory of this 6 // distribution or <https://octave.org/copyright/>. 7 // 8 // This file is part of Octave. 9 // 10 // Octave is free software: you can redistribute it and/or modify it 11 // under the terms of the GNU General Public License as published by 12 // the Free Software Foundation, either version 3 of the License, or 13 // (at your option) any later version. 14 // 15 // Octave is distributed in the hope that it will be useful, but 16 // WITHOUT ANY WARRANTY; without even the implied warranty of 17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18 // GNU General Public License for more details. 19 // 20 // You should have received a copy of the GNU General Public License 21 // along with Octave; see the file COPYING. If not, see 22 // <https://www.gnu.org/licenses/>. 23 // 24 //////////////////////////////////////////////////////////////////////// 25 26 #if defined (HAVE_CONFIG_H) 27 # include "config.h" 28 #endif 29 30 #include "defun.h" 31 #include "error.h" 32 #include "errwarn.h" 33 #include "ovl.h" 34 #include "ops.h" 35 #include "ov-re-diag.h" 36 #include "ov-cx-diag.h" 37 #include "ov-flt-re-diag.h" 38 #include "ov-flt-cx-diag.h" 39 #include "ov-perm.h" 40 41 DEFUN (pinv, args, , 42 doc: /* -*- texinfo -*- 43 @deftypefn {} {} pinv (@var{x}) 44 @deftypefnx {} {} pinv (@var{x}, @var{tol}) 45 Return the @nospell{Moore-Penrose} pseudoinverse of @var{x}. 46 47 Singular values less than @var{tol} are ignored. 48 49 If the second argument is omitted, it is taken to be 50 51 @example 52 tol = max ([rows(@var{x}), columns(@var{x})]) * norm (@var{x}) * eps 53 @end example 54 55 @seealso{inv, ldivide} 56 @end deftypefn */) 57 { 58 int nargin = args.length (); 59 60 if (nargin < 1 || nargin > 2) 61 print_usage (); 62 63 octave_value arg = args(0); 64 65 if (arg.isempty ()) 66 return ovl (Matrix ()); 67 68 octave_value retval; 69 70 bool isfloat = arg.is_single_type (); 71 72 if (arg.is_diag_matrix ()) 73 { 74 if (isfloat) 75 { 76 float tol = 0.0; 77 if (nargin == 2) 78 tol = args(1).float_value (); 79 80 if (tol < 0.0) 81 error ("pinv: TOL must be greater than zero"); 82 83 if (arg.isreal ()) 84 retval = arg.float_diag_matrix_value ().pseudo_inverse (tol); 85 else 86 retval = arg.float_complex_diag_matrix_value ().pseudo_inverse (tol); 87 } 88 else 89 { 90 double tol = 0.0; 91 if (nargin == 2) 92 tol = args(1).double_value (); 93 94 if (tol < 0.0) 95 error ("pinv: TOL must be greater than zero"); 96 97 if (arg.isreal ()) 98 retval = arg.diag_matrix_value ().pseudo_inverse (tol); 99 else 100 retval = arg.complex_diag_matrix_value ().pseudo_inverse (tol); 101 } 102 } 103 else if (arg.is_perm_matrix ()) 104 { 105 retval = arg.perm_matrix_value ().inverse (); 106 } 107 else if (isfloat) 108 { 109 float tol = 0.0; 110 if (nargin == 2) 111 tol = args(1).float_value (); 112 113 if (tol < 0.0) 114 error ("pinv: TOL must be greater than zero"); 115 116 if (arg.isreal ()) 117 { 118 FloatMatrix m = arg.float_matrix_value (); 119 120 retval = m.pseudo_inverse (tol); 121 } 122 else if (arg.iscomplex ()) 123 { 124 FloatComplexMatrix m = arg.float_complex_matrix_value (); 125 126 retval = m.pseudo_inverse (tol); 127 } 128 else 129 err_wrong_type_arg ("pinv", arg); 130 } 131 else 132 { 133 double tol = 0.0; 134 if (nargin == 2) 135 tol = args(1).double_value (); 136 137 if (tol < 0.0) 138 error ("pinv: TOL must be greater than zero"); 139 140 if (arg.isreal ()) 141 { 142 Matrix m = arg.matrix_value (); 143 144 retval = m.pseudo_inverse (tol); 145 } 146 else if (arg.iscomplex ()) 147 { 148 ComplexMatrix m = arg.complex_matrix_value (); 149 150 retval = m.pseudo_inverse (tol); 151 } 152 else 153 err_wrong_type_arg ("pinv", arg); 154 } 155 156 return retval; 157 } 158 159 /* 160 %!shared a, b, tol, hitol, d, u, x, y 161 %! old_state = rand ("state"); 162 %! restore_state = onCleanup (@() rand ("state", old_state)); 163 %! rand ("state", 42); # initialize generator to make behavior reproducible 164 %! a = reshape (rand*[1:16], 4, 4); # Rank 2 matrix 165 %! b = pinv (a); 166 %! tol = 4e-14; 167 %! hitol = 40*sqrt (eps); 168 %! d = diag ([rand, rand, hitol, hitol]); 169 %! u = rand (4); # Could be singular by freak accident 170 %! x = inv (u)*d*u; 171 %! y = pinv (x, sqrt (eps)); 172 173 ## Verify Penrose conditions for pseudoinverse 174 %!assert (a*b*a, a, tol) 175 %!assert (b*a*b, b, tol) 176 %!assert ((b*a)', b*a, tol) 177 %!assert ((a*b)', a*b, tol) 178 %!assert (x*y*x, x, -hitol) 179 %!assert (y*x*y, y, -hitol) 180 %!assert ((x*y)', x*y, hitol) 181 %!assert ((y*x)', y*x, hitol) 182 183 ## Clear shared variables 184 %!shared 185 186 ## Test pinv for Diagonal matrices 187 %!test 188 %! x = diag ([3 2 1 0 -0.5]); 189 %! y = pinv (x); 190 %! assert (typeinfo (y)(1:8), "diagonal"); 191 %! assert (isa (y, "double")); 192 %! assert (diag (y), [1/3, 1/2, 1, 0 1/-0.5]'); 193 %! y = pinv (x, 1); 194 %! assert (diag (y), [1/3 1/2 1 0 0]'); 195 %! y = pinv (x, 2); 196 %! assert (diag (y), [1/3 1/2 0 0 0]'); 197 198 ## Test special case of 0 scalars and vectors 199 %!assert (pinv (0), 0) 200 %!assert (pinv ([0, 0, 0]), [0; 0; 0]) 201 %!assert (pinv (single (0)), single (0)) 202 %!assert (pinv (single ([0, 0, 0])), single ([0; 0; 0])) 203 %!assert (pinv (complex (0,0)), 0) 204 %!assert (pinv (complex ([0,0,0], [0,0,0])), [0; 0; 0]) 205 %!assert (pinv (complex (single (0),0)), single (0)) 206 %!assert (pinv (complex (single ([0,0,0]), [0,0,0])), single ([0; 0; 0])) 207 */ 208