Mercurial > octave
view libinterp/corefcn/qr.cc @ 30241:ba4aebad10d1
Return logical status variable from file functions for Matlab compatibility.
* NEWS: Announce change.
* dirfns.cc (Frmdir): Update documentation. Change return status to true/false
rather than 1.0/0.0.
* copyfile.m: Update documentation. Change internal variable "sts" to logical
variable. Change any assignments to status to true/false.
* movefile.m: Update documentation. Change internal variable "sts" to logical
variable. Change any assignments to status to true/false.
* mkdir.m: Update documentation.
* mkdir.m (mkdir_recur): Change internal variable "status" to logical.
author | Rik <rik@octave.org> |
---|---|
date | Wed, 13 Oct 2021 10:36:20 -0700 |
parents | 7d6709900da7 |
children | 91c6288781ba c980f938fdd8 |
line wrap: on
line source
//////////////////////////////////////////////////////////////////////// // // Copyright (C) 1996-2021 The Octave Project Developers // // See the file COPYRIGHT.md in the top-level directory of this // distribution or <https://octave.org/copyright/>. // // This file is part of Octave. // // Octave is free software: you can redistribute it and/or modify it // under the terms of the GNU General Public License as published by // the Free Software Foundation, either version 3 of the License, or // (at your option) any later version. // // Octave is distributed in the hope that it will be useful, but // WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU General Public License for more details. // // You should have received a copy of the GNU General Public License // along with Octave; see the file COPYING. If not, see // <https://www.gnu.org/licenses/>. // //////////////////////////////////////////////////////////////////////// #if defined (HAVE_CONFIG_H) # include "config.h" #endif #include <string> #include "MArray.h" #include "Matrix.h" #include "qr.h" #include "qrp.h" #include "sparse-qr.h" #include "defun.h" #include "error.h" #include "errwarn.h" #include "ov.h" #include "ovl.h" OCTAVE_NAMESPACE_BEGIN /* ## Restore all rand* "state" values %!function restore_rand_states (state) %! rand ("state", state.rand); %! randn ("state", state.randn); %!endfunction %!shared old_state, restore_state %! ## Save and restore the states of both random number generators that are %! ## tested by the unit tests in this file. %! old_state.rand = rand ("state"); %! old_state.randn = randn ("state"); %! restore_state = onCleanup (@() restore_rand_states (old_state)); */ template <typename MT> static octave_value get_qr_r (const math::qr<MT>& fact) { MT R = fact.R (); if (R.issquare () && fact.regular ()) return octave_value (R, MatrixType (MatrixType::Upper)); else return R; } template <typename T> static typename math::qr<T>::type qr_type (int nargout, bool economy) { if (nargout == 0 || nargout == 1) return math::qr<T>::raw; else if (economy) return math::qr<T>::economy; else return math::qr<T>::std; } // dense X // // [Q, R] = qr (X): form Q unitary and R upper triangular such // that Q * R = X // // [Q, R] = qr (X, 0): form the economy decomposition such that if X is // m by n then only the first n columns of Q are // computed. // // [Q, R, P] = qr (X): form QRP factorization of X where // P is a permutation matrix such that // A * P = Q * R // // [Q, R, P] = qr (X, 0): form the economy decomposition with // permutation vector P such that Q * R = X(:, P) // // qr (X) alone returns the output of the LAPACK routine dgeqrf, such // that R = triu (qr (X)) // // sparse X // // X = qr (A, B): if M < N, X is the minimum 2-norm solution of // A\B. If M >= N, X is the least squares // approximation of A\B. X is calculated by // SPQR-function SuiteSparseQR_min2norm. // DEFUN (qr, args, nargout, doc: /* -*- texinfo -*- @deftypefn {} {[@var{Q}, @var{R}] =} qr (@var{A}) @deftypefnx {} {[@var{Q}, @var{R}, @var{P}] =} qr (@var{A}) @deftypefnx {} {@var{X} =} qr (@var{A}) # non-sparse A @deftypefnx {} {@var{R} =} qr (@var{A}) # sparse A @deftypefnx {} {@var{X} =} qr (@var{A}, @var{B}) # sparse A @deftypefnx {} {[@var{C}, @var{R}] =} qr (@var{A}, @var{B}) @deftypefnx {} {[@dots{}] =} qr (@dots{}, 0) @deftypefnx {} {[@dots{}] =} qr (@dots{}, "vector") @deftypefnx {} {[@dots{}] =} qr (@dots{}, "matrix") @cindex QR factorization Compute the QR@tie{}factorization of @var{A}, using standard @sc{lapack} subroutines. The QR@tie{}factorization is @tex $QR = A$ where $Q$ is an orthogonal matrix and $R$ is upper triangular. @end tex @ifnottex @example @var{Q} * @var{R} = @var{A} @end example @noindent where @var{Q} is an orthogonal matrix and @var{R} is upper triangular. @end ifnottex For example, given the matrix @code{@var{A} = [1, 2; 3, 4]}, @example [@var{Q}, @var{R}] = qr (@var{A}) @end example @noindent returns @example @group @var{Q} = -0.31623 -0.94868 -0.94868 0.31623 @var{R} = -3.16228 -4.42719 0.00000 -0.63246 @end group @end example @noindent which multiplied together return the original matrix @example @group @var{Q} * @var{R} @result{} 1.0000 2.0000 3.0000 4.0000 @end group @end example If just a single return value is requested then it is either @var{R}, if @var{A} is sparse, or @var{X}, such that @code{@var{R} = triu (@var{X})} if @var{A} is full. (Note: unlike most commands, the single return value is not the first return value when multiple values are requested.) If a third output @var{P} is requested, then @code{qr} calculates the permuted QR@tie{}factorization @tex $QR = AP$ where $Q$ is an orthogonal matrix, $R$ is upper triangular, and $P$ is a permutation matrix. @end tex @ifnottex @example @var{Q} * @var{R} = @var{A} * @var{P} @end example @noindent where @var{Q} is an orthogonal matrix, @var{R} is upper triangular, and @var{P} is a permutation matrix. @end ifnottex If @var{A} is dense, the permuted QR@tie{}factorization has the additional property that the diagonal entries of @var{R} are ordered by decreasing magnitude. In other words, @code{abs (diag (@var{R}))} will be ordered from largest to smallest. If @var{A} is sparse, @var{P} is a fill-reducing ordering of the columns of @var{A}. In that case, the diagonal entries of @var{R} are not ordered by decreasing magnitude. For example, given the matrix @code{@var{A} = [1, 2; 3, 4]}, @example [@var{Q}, @var{R}, @var{P}] = qr (@var{A}) @end example @noindent returns @example @group @var{Q} = -0.44721 -0.89443 -0.89443 0.44721 @var{R} = -4.47214 -3.13050 0.00000 0.44721 @var{P} = 0 1 1 0 @end group @end example If the input matrix @var{A} is sparse, the sparse QR@tie{}factorization is computed by using @sc{SPQR} or @sc{CXSparse} (e.g., if @sc{SPQR} is not available). Because the matrix @var{Q} is, in general, a full matrix, it is recommended to request only one return value @var{R}. In that case, the computation avoids the construction of @var{Q} and returns a sparse @var{R} such that @code{@var{R} = chol (@var{A}' * @var{A})}. If @var{A} is dense, an additional matrix @var{B} is supplied and two return values are requested, then @code{qr} returns @var{C}, where @code{@var{C} = @var{Q}' * @var{B}}. This allows the least squares approximation of @code{@var{A} \ @var{B}} to be calculated as @example @group [@var{C}, @var{R}] = qr (@var{A}, @var{B}) @var{X} = @var{R} \ @var{C} @end group @end example If @var{A} is a sparse MxN matrix and an additional matrix @var{B} is supplied, one or two return values are possible. If one return value @var{X} is requested and M < N, then @var{X} is the minimum 2-norm solution of @w{@code{@var{A} \ @var{B}}}. If M >= N, @var{X} is the least squares approximation @w{of @code{@var{A} \ @var{B}}}. If two return values are requested, @var{C} and @var{R} have the same meaning as in the dense case (@var{C} is dense and @var{R} is sparse). The version with one return parameter should be preferred because it uses less memory and can handle rank-deficient matrices better. If the final argument is the string @qcode{"vector"} then @var{P} is a permutation vector (of the columns of @var{A}) instead of a permutation matrix. In this case, the defining relationship is: @example @var{Q} * @var{R} = @var{A}(:, @var{P}) @end example The default, however, is to return a permutation matrix and this may be explicitly specified by using a final argument of @qcode{"matrix"}. If the final argument is the scalar 0 an @qcode{"economy"} factorization is returned. If the original matrix @var{A} has size MxN and M > N, then the @qcode{"economy"} factorization will calculate just N rows in @var{R} and N columns in @var{Q} and omit the zeros in @var{R}. If M @leq{} N, there is no difference between the economy and standard factorizations. When calculating an @qcode{"economy"} factorization and @var{A} is dense, the output @var{P} is always a vector rather than a matrix. If @var{A} is sparse, output @var{P} is a sparse permutation matrix. Background: The QR factorization has applications in the solution of least squares problems @tex $$ \min_x \left\Vert A x - b \right\Vert_2 $$ @end tex @ifnottex @example min norm (A*x - b) @end example @end ifnottex for overdetermined systems of equations (i.e., @tex $A$ @end tex @ifnottex @var{A} @end ifnottex is a tall, thin matrix). The permuted QR@tie{}factorization @code{[@var{Q}, @var{R}, @var{P}] = qr (@var{A})} allows the construction of an orthogonal basis of @code{span (A)}. @seealso{chol, hess, lu, qz, schur, svd, qrupdate, qrinsert, qrdelete, qrshift} @end deftypefn */) { int nargin = args.length (); if (nargin < 1 || nargin > 3) print_usage (); octave_value_list retval; octave_value arg = args(0); bool economy = false; bool is_cmplx = false; bool have_b = 0; bool vector_p = 0; if (arg.iscomplex ()) is_cmplx = true; if (nargin > 1) { have_b = true; if (args(nargin-1).is_scalar_type ()) { int val = args(nargin-1).int_value (); if (val == 0) { economy = true; have_b = (nargin > 2); } else if (nargin == 3) // argument 3 should be 0 or a string print_usage (); } else if (args(nargin-1).is_string ()) { std::string str = args(nargin-1).string_value (); if (str == "vector") vector_p = true; else if (str != "matrix") error ("qr: type for P must be 'matrix' or 'vector', not %s", str.c_str ()); have_b = (nargin > 2); } else if (! args(nargin-1).is_matrix_type ()) err_wrong_type_arg ("qr", args(nargin-1)); else if (nargin == 3) // should be caught by is_scalar_type or is_string print_usage (); if (have_b && args(1).iscomplex ()) is_cmplx = true; } if (arg.issparse ()) { if (nargout > 3) error ("qr: too many output arguments"); if (is_cmplx) { if (have_b && nargout == 1) { octave_idx_type info; if (! args(1).issparse () && args(1).iscomplex ()) retval = ovl (math::sparse_qr<SparseComplexMatrix>::solve <MArray<Complex>, ComplexMatrix> (arg.sparse_complex_matrix_value (), args(1).complex_matrix_value (), info)); else if (args(1).issparse () && args(1).iscomplex ()) retval = ovl (math::sparse_qr<SparseComplexMatrix>::solve <SparseComplexMatrix, SparseComplexMatrix> (arg.sparse_complex_matrix_value (), args(1).sparse_complex_matrix_value (), info)); else if (! args(1).issparse () && ! args(1).iscomplex ()) retval = ovl (math::sparse_qr<SparseComplexMatrix>::solve <MArray<double>, ComplexMatrix> (arg.sparse_complex_matrix_value (), args(1).matrix_value (), info)); else if (args(1).issparse () && ! args(1).iscomplex ()) retval = ovl (math::sparse_qr<SparseComplexMatrix>::solve <SparseMatrix, SparseComplexMatrix> (arg.sparse_complex_matrix_value (), args(1).sparse_matrix_value (), info)); else error ("qr: b is not valid"); } else if (have_b && nargout == 2) { math::sparse_qr<SparseComplexMatrix> q (arg.sparse_complex_matrix_value (), 0); retval = ovl (q.C (args(1).complex_matrix_value (), economy), q.R (economy)); } else if (have_b && nargout == 3) { math::sparse_qr<SparseComplexMatrix> q (arg.sparse_complex_matrix_value ()); if (vector_p) retval = ovl (q.C (args(1).complex_matrix_value (), economy), q.R (economy), q.E ()); else retval = ovl (q.C (args(1).complex_matrix_value (), economy), q.R (economy), q.E_MAT ()); } else { if (nargout > 2) { math::sparse_qr<SparseComplexMatrix> q (arg.sparse_complex_matrix_value ()); if (vector_p) retval = ovl (q.Q (economy), q.R (economy), q.E ()); else retval = ovl (q.Q (economy), q.R (economy), q.E_MAT ()); } else if (nargout > 1) { math::sparse_qr<SparseComplexMatrix> q (arg.sparse_complex_matrix_value (), 0); retval = ovl (q.Q (economy), q.R (economy)); } else { math::sparse_qr<SparseComplexMatrix> q (arg.sparse_complex_matrix_value (), 0); retval = ovl (q.R (economy)); } } } else { if (have_b && nargout == 1) { octave_idx_type info; if (args(1).issparse () && ! args(1).iscomplex ()) retval = ovl (math::sparse_qr<SparseMatrix>::solve <SparseMatrix, SparseMatrix> (arg.sparse_matrix_value (), args (1).sparse_matrix_value (), info)); else if (! args(1).issparse () && args(1).iscomplex ()) retval = ovl (math::sparse_qr<SparseMatrix>::solve <MArray<Complex>, ComplexMatrix> (arg.sparse_matrix_value (), args (1).complex_matrix_value (), info)); else if (! args(1).issparse () && ! args(1).iscomplex ()) retval = ovl (math::sparse_qr<SparseMatrix>::solve <MArray<double>, Matrix> (arg.sparse_matrix_value (), args (1).matrix_value (), info)); else if (args(1).issparse () && args(1).iscomplex ()) retval = ovl (math::sparse_qr<SparseMatrix>::solve <SparseComplexMatrix, SparseComplexMatrix> (arg.sparse_matrix_value (), args(1).sparse_complex_matrix_value (), info)); else error ("qr: b is not valid"); } else if (have_b && nargout == 2) { math::sparse_qr<SparseMatrix> q (arg.sparse_matrix_value (), 0); retval = ovl (q.C (args(1).matrix_value (), economy), q.R (economy)); } else if (have_b && nargout == 3) { math::sparse_qr<SparseMatrix> q (arg.sparse_matrix_value ()); if (vector_p) retval = ovl (q.C (args(1).matrix_value (), economy), q.R (economy), q.E ()); else retval = ovl (q.C (args(1).matrix_value (), economy), q.R (economy), q.E_MAT ()); } else { if (nargout > 2) { math::sparse_qr<SparseMatrix> q (arg.sparse_matrix_value ()); if (vector_p) retval = ovl (q.Q (economy), q.R (economy), q.E ()); else retval = ovl (q.Q (economy), q.R (economy), q.E_MAT ()); } else if (nargout > 1) { math::sparse_qr<SparseMatrix> q (arg.sparse_matrix_value (), 0); retval = ovl (q.Q (economy), q.R (economy)); } else { math::sparse_qr<SparseMatrix> q (arg.sparse_matrix_value (), 0); retval = ovl (q.R (economy)); } } } } else { if (have_b && nargout > 2) error ("qr: too many output arguments for dense A with B"); if (arg.is_single_type ()) { if (arg.isreal ()) { math::qr<FloatMatrix>::type type = qr_type<FloatMatrix> (nargout, economy); FloatMatrix m = arg.float_matrix_value (); switch (nargout) { case 0: case 1: { math::qr<FloatMatrix> fact (m, type); retval = ovl (fact.R ()); } break; case 2: { math::qr<FloatMatrix> fact (m, type); retval = ovl (fact.Q (), get_qr_r (fact)); if (have_b) { if (is_cmplx) retval(0) = fact.Q ().transpose () * args(1).float_complex_matrix_value (); else retval(0) = fact.Q ().transpose () * args(1).float_matrix_value (); } } break; default: { math::qrp<FloatMatrix> fact (m, type); if (economy || vector_p) retval = ovl (fact.Q (), get_qr_r (fact), fact.Pvec ()); else retval = ovl (fact.Q (), get_qr_r (fact), fact.P ()); } break; } } else if (arg.iscomplex ()) { math::qr<FloatComplexMatrix>::type type = qr_type<FloatComplexMatrix> (nargout, economy); FloatComplexMatrix m = arg.float_complex_matrix_value (); switch (nargout) { case 0: case 1: { math::qr<FloatComplexMatrix> fact (m, type); retval = ovl (fact.R ()); } break; case 2: { math::qr<FloatComplexMatrix> fact (m, type); retval = ovl (fact.Q (), get_qr_r (fact)); if (have_b) retval(0) = conj (fact.Q ().transpose ()) * args(1).float_complex_matrix_value (); } break; default: { math::qrp<FloatComplexMatrix> fact (m, type); if (economy || vector_p) retval = ovl (fact.Q (), get_qr_r (fact), fact.Pvec ()); else retval = ovl (fact.Q (), get_qr_r (fact), fact.P ()); } break; } } } else { if (arg.isreal ()) { math::qr<Matrix>::type type = qr_type<Matrix> (nargout, economy); Matrix m = arg.matrix_value (); switch (nargout) { case 0: case 1: { math::qr<Matrix> fact (m, type); retval = ovl (fact.R ()); } break; case 2: { math::qr<Matrix> fact (m, type); retval = ovl (fact.Q (), get_qr_r (fact)); if (have_b) { if (is_cmplx) retval(0) = fact.Q ().transpose () * args(1).complex_matrix_value (); else retval(0) = fact.Q ().transpose () * args(1).matrix_value (); } } break; default: { math::qrp<Matrix> fact (m, type); if (economy || vector_p) retval = ovl (fact.Q (), get_qr_r (fact), fact.Pvec ()); else retval = ovl (fact.Q (), get_qr_r (fact), fact.P ()); } break; } } else if (arg.iscomplex ()) { math::qr<ComplexMatrix>::type type = qr_type<ComplexMatrix> (nargout, economy); ComplexMatrix m = arg.complex_matrix_value (); switch (nargout) { case 0: case 1: { math::qr<ComplexMatrix> fact (m, type); retval = ovl (fact.R ()); } break; case 2: { math::qr<ComplexMatrix> fact (m, type); retval = ovl (fact.Q (), get_qr_r (fact)); if (have_b) retval(0) = conj (fact.Q ().transpose ()) * args(1).complex_matrix_value (); } break; default: { math::qrp<ComplexMatrix> fact (m, type); if (economy || vector_p) retval = ovl (fact.Q (), get_qr_r (fact), fact.Pvec ()); else retval = ovl (fact.Q (), get_qr_r (fact), fact.P ()); } break; } } else err_wrong_type_arg ("qr", arg); } } return retval; } /* %!test %! a = [0, 2, 1; 2, 1, 2]; %! %! [q, r] = qr (a); %! [qe, re] = qr (a, 0); %! %! assert (q * r, a, sqrt (eps)); %! assert (qe * re, a, sqrt (eps)); %!test %! a = [0, 2, 1; 2, 1, 2]; %! %! [q, r] = qr (a); %! [qe, re] = qr (a, 0); %! %! assert (q * r, a, sqrt (eps)); %! assert (qe * re, a, sqrt (eps)); %!test %! a = [0, 2, 1; 2, 1, 2]; %! %! [q, r, p] = qr (a); # FIXME: not giving right dimensions. %! [qe, re, pe] = qr (a, 0); %! %! assert (q * r, a * p, sqrt (eps)); %! assert (qe * re, a(:, pe), sqrt (eps)); %!test %! a = [0, 2; 2, 1; 1, 2]; %! %! [q, r] = qr (a); %! [qe, re] = qr (a, 0); %! %! assert (q * r, a, sqrt (eps)); %! assert (qe * re, a, sqrt (eps)); %!test %! a = [0, 2; 2, 1; 1, 2]; %! %! [q, r, p] = qr (a); %! [qe, re, pe] = qr (a, 0); %! %! assert (q * r, a * p, sqrt (eps)); %! assert (qe * re, a(:, pe), sqrt (eps)); %!test %! a = [0, 2, 1; 2, 1, 2; 3, 1, 2]; %! b = [1, 3, 2; 1, 1, 0; 3, 0, 2]; %! %! [q, r] = qr (a); %! [c, re] = qr (a, b); %! %! assert (r, re, sqrt (eps)); %! assert (q'*b, c, sqrt (eps)); %!test %! a = [0, 2, i; 2, 1, 2; 3, 1, 2]; %! b = [1, 3, 2; 1, i, 0; 3, 0, 2]; %! %! [q, r] = qr (a); %! [c, re] = qr (a, b); %! %! assert (r, re, sqrt (eps)); %! assert (q'*b, c, sqrt (eps)); %!test %! a = [0, 2, i; 2, 1, 2; 3, 1, 2]; %! b = [1, 3, 2; 1, 1, 0; 3, 0, 2]; %! %! [q, r] = qr (a); %! [c, re] = qr (a, b); %! %! assert (r, re, sqrt (eps)); %! assert (q'*b, c, sqrt (eps)); %!test %! a = [0, 2, 1; 2, 1, 2; 3, 1, 2]; %! b = [1, 3, 2; 1, i, 0; 3, 0, 2]; %! %! [q, r] = qr (a); %! [c, re] = qr (a, b); %! %! assert (r, re, sqrt (eps)); %! assert (q'*b, c, sqrt (eps)); %!test %! assert (qr (zeros (0, 0)), zeros (0, 0)) %! assert (qr (zeros (1, 0)), zeros (1, 0)) %! assert (qr (zeros (0, 1)), zeros (0, 1)) %!error qr () %!error qr ([1, 2; 3, 4], 0, 2) %!error <too many output arguments for dense A with B> %! [q, r, p] = qr (rand (3, 2), rand (3, 1)); %!error <too many output arguments for dense A with B> %! [q, r, p] = qr (rand (3, 2), rand (3, 1), 0); %!function retval = __testqr (q, r, a, p) %! tol = 100*eps (class (q)); %! retval = 0; %! if (nargin == 3) %! n1 = norm (q*r - a); %! n2 = norm (q'*q - eye (columns (q))); %! retval = (n1 < tol && n2 < tol); %! else %! n1 = norm (q'*q - eye (columns (q))); %! retval = (n1 < tol); %! if (isvector (p)) %! n2 = norm (q*r - a(:,p)); %! retval = (retval && n2 < tol); %! else %! n2 = norm (q*r - a*p); %! retval = (retval && n2 < tol); %! endif %! endif %!endfunction %!test %! t = ones (24, 1); %! j = 1; %! %! if (false) # eliminate big matrix tests %! a = rand (5000, 20); %! [q,r] = qr (a, 0); t(j++) = __testqr (q, r, a); %! [q,r] = qr (a',0); t(j++) = __testqr (q, r, a'); %! [q,r,p] = qr (a, 0); t(j++) = __testqr (q, r, a, p); %! [q,r,p] = qr (a',0); t(j++) = __testqr (q, r, a', p); %! %! a = a+1i*eps; %! [q,r] = qr (a, 0); t(j++) = __testqr (q, r, a); %! [q,r] = qr (a',0); t(j++) = __testqr (q, r, a'); %! [q,r,p] = qr (a, 0); t(j++) = __testqr (q, r, a, p); %! [q,r,p] = qr (a',0); t(j++) = __testqr (q, r, a', p); %! endif %! %! a = [ ones(1,15); sqrt(eps)*eye(15) ]; %! [q,r] = qr (a); t(j++) = __testqr (q, r, a); %! [q,r] = qr (a'); t(j++) = __testqr (q, r, a'); %! [q,r,p] = qr (a); t(j++) = __testqr (q, r, a, p); %! [q,r,p] = qr (a'); t(j++) = __testqr (q, r, a', p); %! %! a = a+1i*eps; %! [q,r] = qr (a); t(j++) = __testqr (q, r, a); %! [q,r] = qr (a'); t(j++) = __testqr (q, r, a'); %! [q,r,p] = qr (a); t(j++) = __testqr (q, r, a, p); %! [q,r,p] = qr (a'); t(j++) = __testqr (q, r, a', p); %! %! a = [ ones(1,15); sqrt(eps)*eye(15) ]; %! [q,r] = qr (a, 0); t(j++) = __testqr (q, r, a); %! [q,r] = qr (a',0); t(j++) = __testqr (q, r, a'); %! [q,r,p] = qr (a, 0); t(j++) = __testqr (q, r, a, p); %! [q,r,p] = qr (a',0); t(j++) = __testqr (q, r, a', p); %! %! a = a+1i*eps; %! [q,r] = qr (a, 0); t(j++) = __testqr (q, r, a); %! [q,r] = qr (a',0); t(j++) = __testqr (q, r, a'); %! [q,r,p] = qr (a, 0); t(j++) = __testqr (q, r, a, p); %! [q,r,p] = qr (a',0); t(j++) = __testqr (q, r, a', p); %! %! a = [ 611 196 -192 407 -8 -52 -49 29 %! 196 899 113 -192 -71 -43 -8 -44 %! -192 113 899 196 61 49 8 52 %! 407 -192 196 611 8 44 59 -23 %! -8 -71 61 8 411 -599 208 208 %! -52 -43 49 44 -599 411 208 208 %! -49 -8 8 59 208 208 99 -911 %! 29 -44 52 -23 208 208 -911 99 ]; %! [q,r] = qr (a); %! %! assert (all (t) && norm (q*r - a) < 5000*eps); %!test %! a = single ([0, 2, 1; 2, 1, 2]); %! %! [q, r] = qr (a); %! [qe, re] = qr (a, 0); %! %! assert (q * r, a, sqrt (eps ("single"))); %! assert (qe * re, a, sqrt (eps ("single"))); %!test %! a = single ([0, 2, 1; 2, 1, 2]); %! %! [q, r, p] = qr (a); # FIXME: not giving right dimensions. %! [qe, re, pe] = qr (a, 0); %! %! assert (q * r, a * p, sqrt (eps ("single"))); %! assert (qe * re, a(:, pe), sqrt (eps ("single"))); %!test %! a = single ([0, 2; 2, 1; 1, 2]); %! %! [q, r] = qr (a); %! [qe, re] = qr (a, 0); %! %! assert (q * r, a, sqrt (eps ("single"))); %! assert (qe * re, a, sqrt (eps ("single"))); %!test %! a = single ([0, 2; 2, 1; 1, 2]); %! %! [q, r, p] = qr (a); %! [qe, re, pe] = qr (a, 0); %! %! assert (q * r, a * p, sqrt (eps ("single"))); %! assert (qe * re, a(:, pe), sqrt (eps ("single"))); %!test %! a = single ([0, 2, 1; 2, 1, 2; 3, 1, 2]); %! b = single ([1, 3, 2; 1, 1, 0; 3, 0, 2]); %! %! [q, r] = qr (a); %! [c, re] = qr (a, b); %! %! assert (r, re, sqrt (eps ("single"))); %! assert (q'*b, c, sqrt (eps ("single"))); %!test %! a = single ([0, 2, i; 2, 1, 2; 3, 1, 2]); %! b = single ([1, 3, 2; 1, i, 0; 3, 0, 2]); %! %! [q, r] = qr (a); %! [c, re] = qr (a, b); %! %! assert (r, re, sqrt (eps ("single"))); %! assert (q'*b, c, sqrt (eps ("single"))); %!test %! a = single ([0, 2, i; 2, 1, 2; 3, 1, 2]); %! b = single ([1, 3, 2; 1, 1, 0; 3, 0, 2]); %! %! [q, r] = qr (a); %! [c, re] = qr (a, b); %! %! assert (r, re, sqrt (eps)); %! assert (q'*b, c, sqrt (eps)); %!test %! a = single ([0, 2, 1; 2, 1, 2; 3, 1, 2]); %! b = single ([1, 3, 2; 1, i, 0; 3, 0, 2]); %! %! [q, r] = qr (a); %! [c, re] = qr (a, b); %! %! assert (r, re, sqrt (eps ("single"))); %! assert (q'*b, c, sqrt (eps ("single"))); %!test %! t = ones (24, 1); %! j = 1; %! %! if (false) # eliminate big matrix tests %! a = rand (5000,20); %! [q,r] = qr (a, 0); t(j++) = __testqr (q, r, a); %! [q,r] = qr (a',0); t(j++) = __testqr (q, r, a'); %! [q,r,p] = qr (a, 0); t(j++) = __testqr (q, r, a, p); %! [q,r,p] = qr (a',0); t(j++) = __testqr (q, r, a', p); %! %! a = a+1i*eps ("single"); %! [q,r] = qr (a, 0); t(j++) = __testqr (q, r, a); %! [q,r] = qr (a',0); t(j++) = __testqr (q, r, a'); %! [q,r,p] = qr (a, 0); t(j++) = __testqr (q, r, a, p); %! [q,r,p] = qr (a',0); t(j++) = __testqr (q, r, a', p); %! endif %! %! a = [ ones(1,15); sqrt(eps("single"))*eye(15) ]; %! [q,r] = qr (a); t(j++) = __testqr (q, r, a); %! [q,r] = qr (a'); t(j++) = __testqr (q, r, a'); %! [q,r,p] = qr (a); t(j++) = __testqr (q, r, a, p); %! [q,r,p] = qr (a'); t(j++) = __testqr (q, r, a', p); %! %! a = a+1i*eps ("single"); %! [q,r] = qr (a); t(j++) = __testqr (q, r, a); %! [q,r] = qr (a'); t(j++) = __testqr (q, r, a'); %! [q,r,p] = qr (a); t(j++) = __testqr (q, r, a, p); %! [q,r,p] = qr (a'); t(j++) = __testqr (q, r, a', p); %! %! a = [ ones(1,15); sqrt(eps("single"))*eye(15) ]; %! [q,r] = qr (a, 0); t(j++) = __testqr (q, r, a); %! [q,r] = qr (a',0); t(j++) = __testqr (q, r, a'); %! [q,r,p] = qr (a, 0); t(j++) = __testqr (q, r, a, p); %! [q,r,p] = qr (a',0); t(j++) = __testqr (q, r, a', p); %! %! a = a+1i*eps ("single"); %! [q,r] = qr (a, 0); t(j++) = __testqr (q, r, a); %! [q,r] = qr (a',0); t(j++) = __testqr (q, r, a'); %! [q,r,p] = qr (a, 0); t(j++) = __testqr (q, r, a, p); %! [q,r,p] = qr (a',0); t(j++) = __testqr (q, r, a',p); %! %! a = [ 611 196 -192 407 -8 -52 -49 29 %! 196 899 113 -192 -71 -43 -8 -44 %! -192 113 899 196 61 49 8 52 %! 407 -192 196 611 8 44 59 -23 %! -8 -71 61 8 411 -599 208 208 %! -52 -43 49 44 -599 411 208 208 %! -49 -8 8 59 208 208 99 -911 %! 29 -44 52 -23 208 208 -911 99 ]; %! [q,r] = qr (a); %! %! assert (all (t) && norm (q*r-a) < 5000*eps ("single")); ## The deactivated tests below can't be tested till rectangular back-subs is ## implemented for sparse matrices. %!testif ; (__have_feature__ ("SPQR") && __have_feature__ ("CHOLMOD")) || __have_feature__ ("CXSPARSE") %! n = 20; d = 0.2; %! ## initialize generators to make behavior reproducible %! rand ("state", 42); %! randn ("state", 42); %! a = sprandn (n,n,d) + speye (n,n); %! r = qr (a); %! assert (r'*r, a'*a, 1e-10); %!testif HAVE_COLAMD; (__have_feature__ ("SPQR") && __have_feature__ ("CHOLMOD")) || __have_feature__ ("CXSPARSE") %! n = 20; d = 0.2; %! ## initialize generators to make behavior reproducible %! rand ("state", 42); %! randn ("state", 42); %! a = sprandn (n,n,d) + speye (n,n); %! q = symamd (a); %! a = a(q,q); %! r = qr (a); %! assert (r'*r, a'*a, 1e-10); %!testif ; (__have_feature__ ("SPQR") && __have_feature__ ("CHOLMOD")) || __have_feature__ ("CXSPARSE") %! n = 20; d = 0.2; %! ## initialize generators to make behavior reproducible %! rand ("state", 42); %! randn ("state", 42); %! a = sprandn (n,n,d) + speye (n,n); %! [c,r] = qr (a, ones (n,1)); %! assert (r\c, full (a)\ones (n,1), 10e-10); %!testif ; (__have_feature__ ("SPQR") && __have_feature__ ("CHOLMOD")) || __have_feature__ ("CXSPARSE") %! n = 20; d = 0.2; %! ## initialize generators to make behavior reproducible %! rand ("state", 42); %! randn ("state", 42); %! a = sprandn (n,n,d) + speye (n,n); %! b = randn (n,2); %! [c,r] = qr (a, b); %! assert (r\c, full (a)\b, 10e-10); ## Test under-determined systems!! %!#testif ; (__have_feature__ ("SPQR") && __have_feature__ ("CHOLMOD")) || __have_feature__ ("CXSPARSE") %! n = 20; d = 0.2; %! ## initialize generators to make behavior reproducible %! rand ("state", 42); %! randn ("state", 42); %! a = sprandn (n,n+1,d) + speye (n,n+1); %! b = randn (n,2); %! [c,r] = qr (a, b); %! assert (r\c, full (a)\b, 10e-10); %!testif ; (__have_feature__ ("SPQR") && __have_feature__ ("CHOLMOD")) || __have_feature__ ("CXSPARSE") %! n = 20; d = 0.2; %! ## initialize generators to make behavior reproducible %! rand ("state", 42); %! randn ("state", 42); %! a = 1i*sprandn (n,n,d) + speye (n,n); %! r = qr (a); %! assert (r'*r,a'*a,1e-10); %!testif HAVE_COLAMD; (__have_feature__ ("SPQR") && __have_feature__ ("CHOLMOD")) || __have_feature__ ("CXSPARSE") %! n = 20; d = 0.2; %! ## initialize generators to make behavior reproducible %! rand ("state", 42); %! randn ("state", 42); %! a = 1i*sprandn (n,n,d) + speye (n,n); %! q = symamd (a); %! a = a(q,q); %! r = qr (a); %! assert (r'*r, a'*a, 1e-10); %!testif ; (__have_feature__ ("SPQR") && __have_feature__ ("CHOLMOD")) || __have_feature__ ("CXSPARSE") %! n = 20; d = 0.2; %! ## initialize generators to make behavior reproducible %! rand ("state", 42); %! randn ("state", 42); %! a = 1i*sprandn (n,n,d) + speye (n,n); %! [c,r] = qr (a, ones (n,1)); %! assert (r\c, full (a)\ones (n,1), 10e-10); %!testif ; (__have_feature__ ("SPQR") && __have_feature__ ("CHOLMOD")) || __have_feature__ ("CXSPARSE") %! n = 20; d = 0.2; %! ## initialize generators to make behavior reproducible %! rand ("state", 42); %! randn ("state", 42); %! a = 1i*sprandn (n,n,d) + speye (n,n); %! b = randn (n,2); %! [c,r] = qr (a, b); %! assert (r\c, full (a)\b, 10e-10); ## Test under-determined systems!! %!#testif ; (__have_feature__ ("SPQR") && __have_feature__ ("CHOLMOD")) || __have_feature__ ("CXSPARSE") %! n = 20; d = 0.2; %! ## initialize generators to make behavior reproducible %! rand ("state", 42); %! randn ("state", 42); %! a = 1i*sprandn (n,n+1,d) + speye (n,n+1); %! b = randn (n, 2); %! [c, r] = qr (a, b); %! assert (r\c, full (a)\b, 10e-10); %!testif HAVE_SPQR, HAVE_CHOLMOD %! n = 12; m = 20; d = 0.2; %! ## initialize generators to make behavior reproducible %! rand ("state", 42); %! randn ("state", 42); %! a = sprandn (m, n, d); %! b = randn (m, 2); %! [c, r] = qr (a, b); %! assert (r\c, full (a)\b, 10e-10); %!testif HAVE_SPQR, HAVE_CHOLMOD %! n = 12; m = 20; d = 0.2; %! ## initialize generators to make behavior reproducible %! rand ("state", 42); %! randn ("state", 42); %! a = sprandn (m, n, d); %! b = sprandn (m, 2, d); %! [c, r] = qr (a, b, 0); %! [c2, r2] = qr (full (a), full (b), 0); %! assert (r\c, r2\c2, 10e-10); %!testif HAVE_SPQR, HAVE_CHOLMOD %! n = 12; m = 20; d = 0.2; %! ## initialize generators to make behavior reproducible %! rand ("state", 42); %! randn ("state", 42); %! a = sprandn (m, n, d); %! b = randn (m, 2); %! [c, r, p] = qr (a, b, "matrix"); %! x = p * (r\c); %! [c2, r2] = qr (full (a), b); %! x2 = r2 \ c2; %! assert (x, x2, 10e-10); %!testif HAVE_SPQR, HAVE_CHOLMOD %! n = 12; m = 20; d = 0.2; %! ## initialize generators to make behavior reproducible %! rand ("state", 42); %! randn ("state", 42); %! a = sprandn (m, n, d); %! [q, r, p] = qr (a, "matrix"); %! assert (q * r, a * p, 10e-10); %!testif HAVE_SPQR, HAVE_CHOLMOD %! n = 12; m = 20; d = 0.2; %! ## initialize generators to make behavior reproducible %! rand ("state", 42); %! randn ("state", 42); %! a = sprandn (m, n, d); %! b = randn (m, 2); %! x = qr (a, b); %! [c2, r2] = qr (full (a), b); %! assert (x, r2\c2, 10e-10); %!testif HAVE_SPQR, HAVE_CHOLMOD %! n = 12; m = 20; d = 0.2; %! ## initialize generators to make behavior reproducible %! rand ("state", 42); %! randn ("state", 42); %! a = sprandn (m, n, d); %! b = i * randn (m, 2); %! x = qr (a, b); %! [c2, r2] = qr (full (a), b); %! assert (x, r2\c2, 10e-10); %!#testif HAVE_SPQR, HAVE_CHOLMOD %! n = 12; m = 20; d = 0.2; %! ## initialize generators to make behavior reproducible %! rand ("state", 42); %! randn ("state", 42); %! a = sprandn (m, n, d); %! b = i * randn (m, 2); %! [c, r] = qr (a, b); %! [c2, r2] = qr (full (a), b); %! assert (r\c, r2\c2, 10e-10); %!testif HAVE_SPQR, HAVE_CHOLMOD %! n = 12; m = 20; d = 0.2; %! ## initialize generators to make behavior reproducible %! rand ("state", 42); %! randn ("state", 42); %! a = sprandn (m, n, d); %! b = i * randn (m, 2); %! [c, r, p] = qr (a, b, "matrix"); %! x = p * (r\c); %! [c2, r2] = qr (full (a), b); %! x2 = r2 \ c2; %! assert (x, x2, 10e-10); %!testif HAVE_SPQR, HAVE_CHOLMOD %! n = 12; m = 20; d = 0.2; %! ## initialize generators to make behavior reproducible %! rand ("state", 42); %! randn ("state", 42); %! a = i * sprandn (m, n, d); %! b = sprandn (m, 2, d); %! [c, r] = qr (a, b, 0); %! [c2, r2] = qr (full (a), full (b), 0); %! assert (r\c, r2\c2, 10e-10); %!testif HAVE_SPQR, HAVE_CHOLMOD %! n = 12; m = 20; d = 0.2; %! ## initialize generators to make behavior reproducible %! rand ("state", 42); %! randn ("state", 42); %! a = i * sprandn (m, n, d); %! b = randn (m, 2); %! [c, r, p] = qr (a, b, "matrix"); %! x = p * (r\c); %! [c2, r2] = qr (full (a), b); %! x2 = r2 \ c2; %! assert(x, x2, 10e-10); %!testif HAVE_SPQR, HAVE_CHOLMOD %! n = 12; m = 20; d = 0.2; %! ## initialize generators to make behavior reproducible %! rand ("state", 42); %! randn ("state", 42); %! a = i * sprandn (m, n, d); %! [q, r, p] = qr (a, "matrix"); %! assert(q * r, a * p, 10e-10); %!testif HAVE_SPQR, HAVE_CHOLMOD %! n = 12; m = 20; d = 0.2; %! ## initialize generators to make behavior reproducible %! rand ("state", 42); %! randn ("state", 42); %! a = i * sprandn (m, n, d); %! b = randn (m, 2); %! x = qr (a, b); %! [c2, r2] = qr (full (a), b); %! assert (x, r2\c2, 10e-10); %!testif HAVE_SPQR, HAVE_CHOLMOD %! a = sparse (5, 6); %! a(3,1) = 0.8; %! a(2,2) = 1.4; %! a(1,6) = -0.5; %! r = qr (a); %! assert (r'*r, a'*a, 10e-10); */ static bool check_qr_dims (const octave_value& q, const octave_value& r, bool allow_ecf = false) { octave_idx_type m = q.rows (); octave_idx_type k = r.rows (); octave_idx_type n = r.columns (); return ((q.ndims () == 2 && r.ndims () == 2 && k == q.columns ()) && (m == k || (allow_ecf && k == n && k < m))); } static bool check_index (const octave_value& i, bool vector_allowed = false) { return ((i.isreal () || i.isinteger ()) && (i.is_scalar_type () || vector_allowed)); } DEFUN (qrupdate, args, , doc: /* -*- texinfo -*- @deftypefn {} {[@var{Q1}, @var{R1}] =} qrupdate (@var{Q}, @var{R}, @var{u}, @var{v}) Update a QR factorization given update vectors or matrices. Given a QR@tie{}factorization of a real or complex matrix @w{@var{A} = @var{Q}*@var{R}}, @var{Q}@tie{}unitary and @var{R}@tie{}upper trapezoidal, return the QR@tie{}factorization of @w{@var{A} + @var{u}*@var{v}'}, where @var{u} and @var{v} are column vectors (rank-1 update) or matrices with equal number of columns (rank-k update). Notice that the latter case is done as a sequence of rank-1 updates; thus, for k large enough, it will be both faster and more accurate to recompute the factorization from scratch. The QR@tie{}factorization supplied may be either full (Q is square) or economized (R is square). @seealso{qr, qrinsert, qrdelete, qrshift} @end deftypefn */) { octave_value_list retval; if (args.length () != 4) print_usage (); octave_value argq = args(0); octave_value argr = args(1); octave_value argu = args(2); octave_value argv = args(3); if (! argq.isnumeric () || ! argr.isnumeric () || ! argu.isnumeric () || ! argv.isnumeric ()) print_usage (); if (! check_qr_dims (argq, argr, true)) error ("qrupdate: Q and R dimensions don't match"); if (argq.isreal () && argr.isreal () && argu.isreal () && argv.isreal ()) { // all real case if (argq.is_single_type () || argr.is_single_type () || argu.is_single_type () || argv.is_single_type ()) { FloatMatrix Q = argq.float_matrix_value (); FloatMatrix R = argr.float_matrix_value (); FloatMatrix u = argu.float_matrix_value (); FloatMatrix v = argv.float_matrix_value (); math::qr<FloatMatrix> fact (Q, R); fact.update (u, v); retval = ovl (fact.Q (), get_qr_r (fact)); } else { Matrix Q = argq.matrix_value (); Matrix R = argr.matrix_value (); Matrix u = argu.matrix_value (); Matrix v = argv.matrix_value (); math::qr<Matrix> fact (Q, R); fact.update (u, v); retval = ovl (fact.Q (), get_qr_r (fact)); } } else { // complex case if (argq.is_single_type () || argr.is_single_type () || argu.is_single_type () || argv.is_single_type ()) { FloatComplexMatrix Q = argq.float_complex_matrix_value (); FloatComplexMatrix R = argr.float_complex_matrix_value (); FloatComplexMatrix u = argu.float_complex_matrix_value (); FloatComplexMatrix v = argv.float_complex_matrix_value (); math::qr<FloatComplexMatrix> fact (Q, R); fact.update (u, v); retval = ovl (fact.Q (), get_qr_r (fact)); } else { ComplexMatrix Q = argq.complex_matrix_value (); ComplexMatrix R = argr.complex_matrix_value (); ComplexMatrix u = argu.complex_matrix_value (); ComplexMatrix v = argv.complex_matrix_value (); math::qr<ComplexMatrix> fact (Q, R); fact.update (u, v); retval = ovl (fact.Q (), get_qr_r (fact)); } } return retval; } /* %!shared A, u, v, Ac, uc, vc %! A = [0.091364 0.613038 0.999083; %! 0.594638 0.425302 0.603537; %! 0.383594 0.291238 0.085574; %! 0.265712 0.268003 0.238409; %! 0.669966 0.743851 0.445057 ]; %! %! u = [0.85082; %! 0.76426; %! 0.42883; %! 0.53010; %! 0.80683 ]; %! %! v = [0.98810; %! 0.24295; %! 0.43167 ]; %! %! Ac = [0.620405 + 0.956953i 0.480013 + 0.048806i 0.402627 + 0.338171i; %! 0.589077 + 0.658457i 0.013205 + 0.279323i 0.229284 + 0.721929i; %! 0.092758 + 0.345687i 0.928679 + 0.241052i 0.764536 + 0.832406i; %! 0.912098 + 0.721024i 0.049018 + 0.269452i 0.730029 + 0.796517i; %! 0.112849 + 0.603871i 0.486352 + 0.142337i 0.355646 + 0.151496i ]; %! %! uc = [0.20351 + 0.05401i; %! 0.13141 + 0.43708i; %! 0.29808 + 0.08789i; %! 0.69821 + 0.38844i; %! 0.74871 + 0.25821i ]; %! %! vc = [0.85839 + 0.29468i; %! 0.20820 + 0.93090i; %! 0.86184 + 0.34689i ]; %! %!test %! [Q,R] = qr (A); %! [Q,R] = qrupdate (Q, R, u, v); %! assert (norm (vec (Q'*Q - eye (5)), Inf) < 1e1*eps); %! assert (norm (vec (triu (R)-R), Inf) == 0); %! assert (norm (vec (Q*R - A - u*v'), Inf) < norm (A)*1e1*eps); %! %!test %! [Q,R] = qr (Ac); %! [Q,R] = qrupdate (Q, R, uc, vc); %! assert (norm (vec (Q'*Q - eye (5)), Inf) < 1e1*eps); %! assert (norm (vec (triu (R)-R), Inf) == 0); %! assert (norm (vec (Q*R - Ac - uc*vc'), Inf) < norm (Ac)*1e1*eps); %!test %! [Q,R] = qr (single (A)); %! [Q,R] = qrupdate (Q, R, single (u), single (v)); %! assert (norm (vec (Q'*Q - eye (5,"single")), Inf) < 1e1*eps ("single")); %! assert (norm (vec (triu (R)-R), Inf) == 0); %! assert (norm (vec (Q*R - single (A) - single (u)*single (v)'), Inf) < norm (single (A))*1e1*eps ("single")); %! %!test %! [Q,R] = qr (single (Ac)); %! [Q,R] = qrupdate (Q, R, single (uc), single (vc)); %! assert (norm (vec (Q'*Q - eye (5,"single")), Inf) < 1e1*eps ("single")); %! assert (norm (vec (triu (R)-R), Inf) == 0); %! assert (norm (vec (Q*R - single (Ac) - single (uc)*single (vc)'), Inf) < norm (single (Ac))*1e1*eps ("single")); */ DEFUN (qrinsert, args, , doc: /* -*- texinfo -*- @deftypefn {} {[@var{Q1}, @var{R1}] =} qrinsert (@var{Q}, @var{R}, @var{j}, @var{x}, @var{orient}) Update a QR factorization given a row or column to insert in the original factored matrix. Given a QR@tie{}factorization of a real or complex matrix @w{@var{A} = @var{Q}*@var{R}}, @var{Q}@tie{}unitary and @var{R}@tie{}upper trapezoidal, return the QR@tie{}factorization of @w{[A(:,1:j-1) x A(:,j:n)]}, where @var{u} is a column vector to be inserted into @var{A} (if @var{orient} is @qcode{"col"}), or the QR@tie{}factorization of @w{[A(1:j-1,:);x;A(:,j:n)]}, where @var{x} is a row vector to be inserted into @var{A} (if @var{orient} is @qcode{"row"}). The default value of @var{orient} is @qcode{"col"}. If @var{orient} is @qcode{"col"}, @var{u} may be a matrix and @var{j} an index vector resulting in the QR@tie{}factorization of a matrix @var{B} such that @w{B(:,@var{j})} gives @var{u} and @w{B(:,@var{j}) = []} gives @var{A}. Notice that the latter case is done as a sequence of k insertions; thus, for k large enough, it will be both faster and more accurate to recompute the factorization from scratch. If @var{orient} is @qcode{"col"}, the QR@tie{}factorization supplied may be either full (Q is square) or economized (R is square). If @var{orient} is @qcode{"row"}, full factorization is needed. @seealso{qr, qrupdate, qrdelete, qrshift} @end deftypefn */) { int nargin = args.length (); if (nargin < 4 || nargin > 5) print_usage (); octave_value argq = args(0); octave_value argr = args(1); octave_value argj = args(2); octave_value argx = args(3); if (! argq.isnumeric () || ! argr.isnumeric () || ! argx.isnumeric () || (nargin > 4 && ! args(4).is_string ())) print_usage (); std::string orient = (nargin < 5) ? "col" : args(4).string_value (); bool col = (orient == "col"); if (! col && orient != "row") error (R"(qrinsert: ORIENT must be "col" or "row")"); if (! check_qr_dims (argq, argr, col) || (! col && argx.rows () != 1)) error ("qrinsert: dimension mismatch"); if (! check_index (argj, col)) error ("qrinsert: invalid index J"); octave_value_list retval; MArray<octave_idx_type> j = argj.octave_idx_type_vector_value (); octave_idx_type one = 1; if (argq.isreal () && argr.isreal () && argx.isreal ()) { // real case if (argq.is_single_type () || argr.is_single_type () || argx.is_single_type ()) { FloatMatrix Q = argq.float_matrix_value (); FloatMatrix R = argr.float_matrix_value (); FloatMatrix x = argx.float_matrix_value (); math::qr<FloatMatrix> fact (Q, R); if (col) fact.insert_col (x, j-one); else fact.insert_row (x.row (0), j(0)-one); retval = ovl (fact.Q (), get_qr_r (fact)); } else { Matrix Q = argq.matrix_value (); Matrix R = argr.matrix_value (); Matrix x = argx.matrix_value (); math::qr<Matrix> fact (Q, R); if (col) fact.insert_col (x, j-one); else fact.insert_row (x.row (0), j(0)-one); retval = ovl (fact.Q (), get_qr_r (fact)); } } else { // complex case if (argq.is_single_type () || argr.is_single_type () || argx.is_single_type ()) { FloatComplexMatrix Q = argq.float_complex_matrix_value (); FloatComplexMatrix R = argr.float_complex_matrix_value (); FloatComplexMatrix x = argx.float_complex_matrix_value (); math::qr<FloatComplexMatrix> fact (Q, R); if (col) fact.insert_col (x, j-one); else fact.insert_row (x.row (0), j(0)-one); retval = ovl (fact.Q (), get_qr_r (fact)); } else { ComplexMatrix Q = argq.complex_matrix_value (); ComplexMatrix R = argr.complex_matrix_value (); ComplexMatrix x = argx.complex_matrix_value (); math::qr<ComplexMatrix> fact (Q, R); if (col) fact.insert_col (x, j-one); else fact.insert_row (x.row (0), j(0)-one); retval = ovl (fact.Q (), get_qr_r (fact)); } } return retval; } /* %!test %! [Q,R] = qr (A); %! [Q,R] = qrinsert (Q, R, 3, u); %! assert (norm (vec (Q'*Q - eye (5)), Inf) < 1e1*eps); %! assert (norm (vec (triu (R) - R), Inf) == 0); %! assert (norm (vec (Q*R - [A(:,1:2) u A(:,3)]), Inf) < norm (A)*1e1*eps); %!test %! [Q,R] = qr (Ac); %! [Q,R] = qrinsert (Q, R, 3, uc); %! assert (norm (vec (Q'*Q - eye (5)), Inf) < 1e1*eps); %! assert (norm (vec (triu (R) - R), Inf) == 0); %! assert (norm (vec (Q*R - [Ac(:,1:2) uc Ac(:,3)]), Inf) < norm (Ac)*1e1*eps); %!test %! x = [0.85082 0.76426 0.42883 ]; %! %! [Q,R] = qr (A); %! [Q,R] = qrinsert (Q, R, 3, x, "row"); %! assert (norm (vec (Q'*Q - eye (6)), Inf) < 1e1*eps); %! assert (norm (vec (triu (R) - R), Inf) == 0); %! assert (norm (vec (Q*R - [A(1:2,:);x;A(3:5,:)]), Inf) < norm (A)*1e1*eps); %!test %! x = [0.20351 + 0.05401i 0.13141 + 0.43708i 0.29808 + 0.08789i ]; %! %! [Q,R] = qr (Ac); %! [Q,R] = qrinsert (Q, R, 3, x, "row"); %! assert (norm (vec (Q'*Q - eye (6)), Inf) < 1e1*eps); %! assert (norm (vec (triu (R) - R), Inf) == 0); %! assert (norm (vec (Q*R - [Ac(1:2,:);x;Ac(3:5,:)]), Inf) < norm (Ac)*1e1*eps); %!test %! [Q,R] = qr (single (A)); %! [Q,R] = qrinsert (Q, R, 3, single (u)); %! assert (norm (vec (Q'*Q - eye (5,"single")), Inf) < 1e1*eps ("single")); %! assert (norm (vec (triu (R) - R), Inf) == 0); %! assert (norm (vec (Q*R - single ([A(:,1:2) u A(:,3)])), Inf) < norm (single (A))*1e1*eps ("single")); %!test %! [Q,R] = qr (single (Ac)); %! [Q,R] = qrinsert (Q, R, 3, single (uc)); %! assert (norm (vec (Q'*Q - eye (5,"single")), Inf) < 1e1*eps ("single")); %! assert (norm (vec (triu (R) - R), Inf) == 0); %! assert (norm (vec (Q*R - single ([Ac(:,1:2) uc Ac(:,3)])), Inf) < norm (single (Ac))*1e1*eps ("single")); %!test %! x = single ([0.85082 0.76426 0.42883 ]); %! %! [Q,R] = qr (single (A)); %! [Q,R] = qrinsert (Q, R, 3, x, "row"); %! assert (norm (vec (Q'*Q - eye (6,"single")), Inf) < 1e1*eps ("single")); %! assert (norm (vec (triu (R) - R), Inf) == 0); %! assert (norm (vec (Q*R - single ([A(1:2,:);x;A(3:5,:)])), Inf) < norm (single (A))*1e1*eps ("single")); %!test %! x = single ([0.20351 + 0.05401i 0.13141 + 0.43708i 0.29808 + 0.08789i ]); %! %! [Q,R] = qr (single (Ac)); %! [Q,R] = qrinsert (Q, R, 3, x, "row"); %! assert (norm (vec (Q'*Q - eye (6,"single")), Inf) < 1e1*eps ("single")); %! assert (norm (vec (triu (R) - R), Inf) == 0); %! assert (norm (vec (Q*R - single ([Ac(1:2,:);x;Ac(3:5,:)])), Inf) < norm (single (Ac))*1e1*eps ("single")); */ DEFUN (qrdelete, args, , doc: /* -*- texinfo -*- @deftypefn {} {[@var{Q1}, @var{R1}] =} qrdelete (@var{Q}, @var{R}, @var{j}, @var{orient}) Update a QR factorization given a row or column to delete from the original factored matrix. Given a QR@tie{}factorization of a real or complex matrix @w{@var{A} = @var{Q}*@var{R}}, @var{Q}@tie{}unitary and @var{R}@tie{}upper trapezoidal, return the QR@tie{}factorization of @w{[A(:,1:j-1), U, A(:,j:n)]}, where @var{u} is a column vector to be inserted into @var{A} (if @var{orient} is @qcode{"col"}), or the QR@tie{}factorization of @w{[A(1:j-1,:);X;A(:,j:n)]}, where @var{x} is a row @var{orient} is @qcode{"row"}). The default value of @var{orient} is @qcode{"col"}. If @var{orient} is @qcode{"col"}, @var{j} may be an index vector resulting in the QR@tie{}factorization of a matrix @var{B} such that @w{A(:,@var{j}) = []} gives @var{B}. Notice that the latter case is done as a sequence of k deletions; thus, for k large enough, it will be both faster and more accurate to recompute the factorization from scratch. If @var{orient} is @qcode{"col"}, the QR@tie{}factorization supplied may be either full (Q is square) or economized (R is square). If @var{orient} is @qcode{"row"}, full factorization is needed. @seealso{qr, qrupdate, qrinsert, qrshift} @end deftypefn */) { int nargin = args.length (); if (nargin < 3 || nargin > 4) print_usage (); octave_value argq = args(0); octave_value argr = args(1); octave_value argj = args(2); if (! argq.isnumeric () || ! argr.isnumeric () || (nargin > 3 && ! args(3).is_string ())) print_usage (); std::string orient = (nargin < 4) ? "col" : args(3).string_value (); bool col = orient == "col"; if (! col && orient != "row") error (R"(qrdelete: ORIENT must be "col" or "row")"); if (! check_qr_dims (argq, argr, col)) error ("qrdelete: dimension mismatch"); MArray<octave_idx_type> j = argj.octave_idx_type_vector_value (); if (! check_index (argj, col)) error ("qrdelete: invalid index J"); octave_value_list retval; octave_idx_type one = 1; if (argq.isreal () && argr.isreal ()) { // real case if (argq.is_single_type () || argr.is_single_type ()) { FloatMatrix Q = argq.float_matrix_value (); FloatMatrix R = argr.float_matrix_value (); math::qr<FloatMatrix> fact (Q, R); if (col) fact.delete_col (j-one); else fact.delete_row (j(0)-one); retval = ovl (fact.Q (), get_qr_r (fact)); } else { Matrix Q = argq.matrix_value (); Matrix R = argr.matrix_value (); math::qr<Matrix> fact (Q, R); if (col) fact.delete_col (j-one); else fact.delete_row (j(0)-one); retval = ovl (fact.Q (), get_qr_r (fact)); } } else { // complex case if (argq.is_single_type () || argr.is_single_type ()) { FloatComplexMatrix Q = argq.float_complex_matrix_value (); FloatComplexMatrix R = argr.float_complex_matrix_value (); math::qr<FloatComplexMatrix> fact (Q, R); if (col) fact.delete_col (j-one); else fact.delete_row (j(0)-one); retval = ovl (fact.Q (), get_qr_r (fact)); } else { ComplexMatrix Q = argq.complex_matrix_value (); ComplexMatrix R = argr.complex_matrix_value (); math::qr<ComplexMatrix> fact (Q, R); if (col) fact.delete_col (j-one); else fact.delete_row (j(0)-one); retval = ovl (fact.Q (), get_qr_r (fact)); } } return retval; } /* %!test %! AA = [0.091364 0.613038 0.027504 0.999083; %! 0.594638 0.425302 0.562834 0.603537; %! 0.383594 0.291238 0.742073 0.085574; %! 0.265712 0.268003 0.783553 0.238409; %! 0.669966 0.743851 0.457255 0.445057 ]; %! %! [Q,R] = qr (AA); %! [Q,R] = qrdelete (Q, R, 3); %! assert (norm (vec (Q'*Q - eye (5)), Inf) < 16*eps); %! assert (norm (vec (triu (R) - R), Inf) == 0); %! assert (norm (vec (Q*R - [AA(:,1:2) AA(:,4)]), Inf) < norm (AA)*1e1*eps); %! %!test %! AA = [0.364554 + 0.993117i 0.669818 + 0.510234i 0.426568 + 0.041337i 0.847051 + 0.233291i; %! 0.049600 + 0.242783i 0.448946 + 0.484022i 0.141155 + 0.074420i 0.446746 + 0.392706i; %! 0.581922 + 0.657416i 0.581460 + 0.030016i 0.219909 + 0.447288i 0.201144 + 0.069132i; %! 0.694986 + 0.000571i 0.682327 + 0.841712i 0.807537 + 0.166086i 0.192767 + 0.358098i; %! 0.945002 + 0.066788i 0.350492 + 0.642638i 0.579629 + 0.048102i 0.600170 + 0.636938i ] * I; %! %! [Q,R] = qr (AA); %! [Q,R] = qrdelete (Q, R, 3); %! assert (norm (vec (Q'*Q - eye (5)), Inf) < 16*eps); %! assert (norm (vec (triu (R) - R), Inf) == 0); %! assert (norm (vec (Q*R - [AA(:,1:2) AA(:,4)]), Inf) < norm (AA)*1e1*eps); %! %!test %! AA = [0.091364 0.613038 0.027504 0.999083; %! 0.594638 0.425302 0.562834 0.603537; %! 0.383594 0.291238 0.742073 0.085574; %! 0.265712 0.268003 0.783553 0.238409; %! 0.669966 0.743851 0.457255 0.445057 ]; %! %! [Q,R] = qr (AA); %! [Q,R] = qrdelete (Q, R, 3, "row"); %! assert (norm (vec (Q'*Q - eye (4)), Inf) < 1e1*eps); %! assert (norm (vec (triu (R) - R), Inf) == 0); %! assert (norm (vec (Q*R - [AA(1:2,:);AA(4:5,:)]), Inf) < norm (AA)*1e1*eps); %! %!test %! AA = [0.364554 + 0.993117i 0.669818 + 0.510234i 0.426568 + 0.041337i 0.847051 + 0.233291i; %! 0.049600 + 0.242783i 0.448946 + 0.484022i 0.141155 + 0.074420i 0.446746 + 0.392706i; %! 0.581922 + 0.657416i 0.581460 + 0.030016i 0.219909 + 0.447288i 0.201144 + 0.069132i; %! 0.694986 + 0.000571i 0.682327 + 0.841712i 0.807537 + 0.166086i 0.192767 + 0.358098i; %! 0.945002 + 0.066788i 0.350492 + 0.642638i 0.579629 + 0.048102i 0.600170 + 0.636938i ] * I; %! %! [Q,R] = qr (AA); %! [Q,R] = qrdelete (Q, R, 3, "row"); %! assert (norm (vec (Q'*Q - eye (4)), Inf) < 1e1*eps); %! assert (norm (vec (triu (R) - R), Inf) == 0); %! assert (norm (vec (Q*R - [AA(1:2,:);AA(4:5,:)]), Inf) < norm (AA)*1e1*eps); %!test %! AA = single ([0.091364 0.613038 0.027504 0.999083; %! 0.594638 0.425302 0.562834 0.603537; %! 0.383594 0.291238 0.742073 0.085574; %! 0.265712 0.268003 0.783553 0.238409; %! 0.669966 0.743851 0.457255 0.445057 ]); %! %! [Q,R] = qr (AA); %! [Q,R] = qrdelete (Q, R, 3); %! assert (norm (vec (Q'*Q - eye (5,"single")), Inf) < 1e1*eps ("single")); %! assert (norm (vec (triu (R) - R), Inf) == 0); %! assert (norm (vec (Q*R - [AA(:,1:2) AA(:,4)]), Inf) < norm (AA)*1e1*eps ("single")); %! %!test %! AA = single ([0.364554 + 0.993117i 0.669818 + 0.510234i 0.426568 + 0.041337i 0.847051 + 0.233291i; %! 0.049600 + 0.242783i 0.448946 + 0.484022i 0.141155 + 0.074420i 0.446746 + 0.392706i; %! 0.581922 + 0.657416i 0.581460 + 0.030016i 0.219909 + 0.447288i 0.201144 + 0.069132i; %! 0.694986 + 0.000571i 0.682327 + 0.841712i 0.807537 + 0.166086i 0.192767 + 0.358098i; %! 0.945002 + 0.066788i 0.350492 + 0.642638i 0.579629 + 0.048102i 0.600170 + 0.636938i ]) * I; %! %! [Q,R] = qr (AA); %! [Q,R] = qrdelete (Q, R, 3); %! assert (norm (vec (Q'*Q - eye (5,"single")), Inf) < 1e1*eps ("single")); %! assert (norm (vec (triu (R) - R), Inf) == 0); %! assert (norm (vec (Q*R - [AA(:,1:2) AA(:,4)]), Inf) < norm (AA)*1e1*eps ("single")); %! %!test %! AA = single ([0.091364 0.613038 0.027504 0.999083; %! 0.594638 0.425302 0.562834 0.603537; %! 0.383594 0.291238 0.742073 0.085574; %! 0.265712 0.268003 0.783553 0.238409; %! 0.669966 0.743851 0.457255 0.445057 ]); %! %! [Q,R] = qr (AA); %! [Q,R] = qrdelete (Q, R, 3, "row"); %! assert (norm (vec (Q'*Q - eye (4,"single")), Inf) < 1.5e1*eps ("single")); %! assert (norm (vec (triu (R) - R), Inf) == 0); %! assert (norm (vec (Q*R - [AA(1:2,:);AA(4:5,:)]), Inf) < norm (AA)*1e1*eps ("single")); %!testif HAVE_QRUPDATE %! ## Same test as above but with more precicision %! AA = single ([0.091364 0.613038 0.027504 0.999083; %! 0.594638 0.425302 0.562834 0.603537; %! 0.383594 0.291238 0.742073 0.085574; %! 0.265712 0.268003 0.783553 0.238409; %! 0.669966 0.743851 0.457255 0.445057 ]); %! %! [Q,R] = qr (AA); %! [Q,R] = qrdelete (Q, R, 3, "row"); %! assert (norm (vec (Q'*Q - eye (4,"single")), Inf) < 1e1*eps ("single")); %! assert (norm (vec (triu (R) - R), Inf) == 0); %! assert (norm (vec (Q*R - [AA(1:2,:);AA(4:5,:)]), Inf) < norm (AA)*1e1*eps ("single")); %! %!test %! AA = single ([0.364554 + 0.993117i 0.669818 + 0.510234i 0.426568 + 0.041337i 0.847051 + 0.233291i; %! 0.049600 + 0.242783i 0.448946 + 0.484022i 0.141155 + 0.074420i 0.446746 + 0.392706i; %! 0.581922 + 0.657416i 0.581460 + 0.030016i 0.219909 + 0.447288i 0.201144 + 0.069132i; %! 0.694986 + 0.000571i 0.682327 + 0.841712i 0.807537 + 0.166086i 0.192767 + 0.358098i; %! 0.945002 + 0.066788i 0.350492 + 0.642638i 0.579629 + 0.048102i 0.600170 + 0.636938i ]) * I; %! %! [Q,R] = qr (AA); %! [Q,R] = qrdelete (Q, R, 3, "row"); %! assert (norm (vec (Q'*Q - eye (4,"single")), Inf) < 1e1*eps ("single")); %! assert (norm (vec (triu (R) - R), Inf) == 0); %! assert (norm (vec (Q*R - [AA(1:2,:);AA(4:5,:)]), Inf) < norm (AA)*1e1*eps ("single")); */ DEFUN (qrshift, args, , doc: /* -*- texinfo -*- @deftypefn {} {[@var{Q1}, @var{R1}] =} qrshift (@var{Q}, @var{R}, @var{i}, @var{j}) Update a QR factorization given a range of columns to shift in the original factored matrix. Given a QR@tie{}factorization of a real or complex matrix @w{@var{A} = @var{Q}*@var{R}}, @var{Q}@tie{}unitary and @var{R}@tie{}upper trapezoidal, return the QR@tie{}factorization of @w{@var{A}(:,p)}, where @w{p} is the permutation @* @code{p = [1:i-1, shift(i:j, 1), j+1:n]} if @w{@var{i} < @var{j}} @* or @* @code{p = [1:j-1, shift(j:i,-1), i+1:n]} if @w{@var{j} < @var{i}}. @* @seealso{qr, qrupdate, qrinsert, qrdelete} @end deftypefn */) { if (args.length () != 4) print_usage (); octave_value argq = args(0); octave_value argr = args(1); octave_value argi = args(2); octave_value argj = args(3); if (! argq.isnumeric () || ! argr.isnumeric ()) print_usage (); if (! check_qr_dims (argq, argr, true)) error ("qrshift: dimensions mismatch"); octave_idx_type i = argi.idx_type_value (); octave_idx_type j = argj.idx_type_value (); if (! check_index (argi) || ! check_index (argj)) error ("qrshift: invalid index I or J"); octave_value_list retval; if (argq.isreal () && argr.isreal ()) { // all real case if (argq.is_single_type () && argr.is_single_type ()) { FloatMatrix Q = argq.float_matrix_value (); FloatMatrix R = argr.float_matrix_value (); math::qr<FloatMatrix> fact (Q, R); fact.shift_cols (i-1, j-1); retval = ovl (fact.Q (), get_qr_r (fact)); } else { Matrix Q = argq.matrix_value (); Matrix R = argr.matrix_value (); math::qr<Matrix> fact (Q, R); fact.shift_cols (i-1, j-1); retval = ovl (fact.Q (), get_qr_r (fact)); } } else { // complex case if (argq.is_single_type () && argr.is_single_type ()) { FloatComplexMatrix Q = argq.float_complex_matrix_value (); FloatComplexMatrix R = argr.float_complex_matrix_value (); math::qr<FloatComplexMatrix> fact (Q, R); fact.shift_cols (i-1, j-1); retval = ovl (fact.Q (), get_qr_r (fact)); } else { ComplexMatrix Q = argq.complex_matrix_value (); ComplexMatrix R = argr.complex_matrix_value (); math::qr<ComplexMatrix> fact (Q, R); fact.shift_cols (i-1, j-1); retval = ovl (fact.Q (), get_qr_r (fact)); } } return retval; } /* %!test %! AA = A.'; %! i = 2; j = 4; p = [1:i-1, shift(i:j,-1), j+1:5]; %! %! [Q,R] = qr (AA); %! [Q,R] = qrshift (Q, R, i, j); %! assert (norm (vec (Q'*Q - eye (3)), Inf) < 1e1*eps); %! assert (norm (vec (triu (R) - R), Inf) == 0); %! assert (norm (vec (Q*R - AA(:,p)), Inf) < norm (AA)*1e1*eps); %! %! j = 2; i = 4; p = [1:j-1, shift(j:i,+1), i+1:5]; %! %! [Q,R] = qr (AA); %! [Q,R] = qrshift (Q, R, i, j); %! assert (norm (vec (Q'*Q - eye (3)), Inf) < 1e1*eps); %! assert (norm (vec (triu (R) - R), Inf) == 0); %! assert (norm (vec (Q*R - AA(:,p)), Inf) < norm (AA)*1e1*eps); %! %!test %! AA = Ac.'; %! i = 2; j = 4; p = [1:i-1, shift(i:j,-1), j+1:5]; %! %! [Q,R] = qr (AA); %! [Q,R] = qrshift (Q, R, i, j); %! assert (norm (vec (Q'*Q - eye (3)), Inf) < 1e1*eps); %! assert (norm (vec (triu (R) - R), Inf) == 0); %! assert (norm (vec (Q*R - AA(:,p)), Inf) < norm (AA)*1e1*eps); %! %! j = 2; i = 4; p = [1:j-1, shift(j:i,+1), i+1:5]; %! %! [Q,R] = qr (AA); %! [Q,R] = qrshift (Q, R, i, j); %! assert (norm (vec (Q'*Q - eye (3)), Inf) < 1e1*eps); %! assert (norm (vec (triu (R) - R), Inf) == 0); %! assert (norm (vec (Q*R - AA(:,p)), Inf) < norm (AA)*1e1*eps); %!test %! AA = single (A).'; %! i = 2; j = 4; p = [1:i-1, shift(i:j,-1), j+1:5]; %! %! [Q,R] = qr (AA); %! [Q,R] = qrshift (Q, R, i, j); %! assert (norm (vec (Q'*Q - eye (3,"single")), Inf) < 1e1*eps ("single")); %! assert (norm (vec (triu (R) - R), Inf) == 0); %! assert (norm (vec (Q*R - AA(:,p)), Inf) < norm (AA)*1e1*eps ("single")); %! %! j = 2; i = 4; p = [1:j-1, shift(j:i,+1), i+1:5]; %! %! [Q,R] = qr (AA); %! [Q,R] = qrshift (Q, R, i, j); %! assert (norm (vec (Q'*Q - eye (3,"single")), Inf) < 1e1*eps ("single")); %! assert (norm (vec (triu (R) - R), Inf) == 0); %! assert (norm (vec (Q*R - AA(:,p)), Inf) < norm (AA)*1e1*eps ("single")); %! %!test %! AA = single (Ac).'; %! i = 2; j = 4; p = [1:i-1, shift(i:j,-1), j+1:5]; %! %! [Q,R] = qr (AA); %! [Q,R] = qrshift (Q, R, i, j); %! assert (norm (vec (Q'*Q - eye (3,"single")), Inf) < 1e1*eps ("single")); %! assert (norm (vec (triu (R) - R), Inf) == 0); %! assert (norm (vec (Q*R - AA(:,p)), Inf) < norm (AA)*1e1*eps ("single")); %! %! j = 2; i = 4; p = [1:j-1, shift(j:i,+1), i+1:5]; %! %! [Q,R] = qr (AA); %! [Q,R] = qrshift (Q, R, i, j); %! assert (norm (vec (Q'*Q - eye (3,"single")), Inf) < 1e1*eps ("single")); %! assert (norm (vec (triu (R) - R), Inf) == 0); %! assert (norm (vec (Q*R - AA(:,p)), Inf) < norm (AA)*1e1*eps ("single")); */ OCTAVE_NAMESPACE_END