Mercurial > octave

changeset 24371:c9d229f1db04

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Correctly return economy or standard QR factorization (bug #52593). * NEWS: Announce change in factorization. * qr.cc (qr_type): Recode function to determine QR type based on nargout and economy variables. * qr.cc (Fqr): Rewrite docstring. Use temporary variable str to hold input argument and simplify input validation. Change all function calls to qr_type() to use new syntax. Change return argument processing to return a permutation vector for either an economy factorization or "vector" argument given.
author Rik <rik@octave.org>
date Wed, 06 Dec 2017 16:51:19 -0800
parents b1d1229d9e83
children 798b56f0b207
files NEWS libinterp/dldfcn/qr.cc
diffstat 2 files changed, 122 insertions(+), 74 deletions(-) [+]
line wrap: on
line diff
--- a/NEWS	Wed Dec 06 14:09:08 2017 -0800
+++ b/NEWS	Wed Dec 06 16:51:19 2017 -0800
@@ -59,9 +59,12 @@
 
       warning ("off", "Octave:quadcc:RelTol-conversion")
 
- ** The properties "BackgroundColor", "EdgeColor", "LineStyle",
-    "LineWidth", and "Margin" are now implemented for text objects using
-    an OpenGL toolkit.
+ ** The qr function now returns a standard factorization unless
+    explicitly instructed to perform an economy factorization by using a
+    final argument of 0.
+
+ ** Text objects now implement the properties "BackgroundColor",
+    "EdgeColor", "LineStyle", "LineWidth", and "Margin".
 
  ** An initial implementation of alpha transparency has been made for
     patch and surface objects.  Printing to svg and pdf is supported.
--- a/libinterp/dldfcn/qr.cc	Wed Dec 06 14:09:08 2017 -0800
+++ b/libinterp/dldfcn/qr.cc	Wed Dec 06 16:51:19 2017 -0800
@@ -53,14 +53,17 @@
 
 template <typename T>
 static typename octave::math::qr<T>::type
-qr_type (int nargin, int nargout)
+qr_type (int nargout, bool economy)
 {
-  return ((nargout == 0 || nargout == 1)
-          ? octave::math::qr<T>::raw
-          : (nargin == 2) ? octave::math::qr<T>::economy : octave::math::qr<T>::std);
+  if (nargout == 0 || nargout == 1)
+    return octave::math::qr<T>::raw;
+  else if (economy)
+    return octave::math::qr<T>::economy;
+  else
+    return octave::math::qr<T>::std;
 }
 
-// [Q, R] = qr (X):      form Q unitary and R upper triangular such
+// [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
@@ -72,7 +75,7 @@
 //                        A * P = Q * R
 //
 // [Q, R, P] = qr (X, 0): form the economy decomposition with
-//                        permutation vector P such that Q * R = X (:, P)
+//                        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))
@@ -80,23 +83,29 @@
 DEFUN_DLD (qr, args, nargout,
            doc: /* -*- texinfo -*-
 @deftypefn  {} {[@var{Q}, @var{R}] =} qr (@var{A})
-@deftypefnx {} {[@var{Q}, @var{R}, @var{P}] =} qr (@var{A}) # non-sparse A
-@deftypefnx {} {@var{X} =} qr (@var{A})
-@deftypefnx {} {@var{R} =} qr (@var{A}) # sparse A
+@deftypefnx {} {[@var{Q}, @var{R}, @var{P}] =} qr (@var{A})  # non-sparse A
+@deftypefnx {} {@var{X} =} qr (@var{A})  # non-sparse A
+@deftypefnx {} {@var{R} =} qr (@var{A})  # 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')
+@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
-@code{@var{Q} * @var{R} = @var{A}} where @var{Q} is an orthogonal matrix and
-@var{R} is upper triangular.
+
+@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]},
@@ -122,41 +131,45 @@
 @end group
 @end example
 
-The @code{qr} factorization has applications in the solution of least
-squares problems
+@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 the matrix @var{A} is full, and a third output @var{P} is requested, then
+@code{qr} calculates the permuted QR@tie{}factorization
 @tex
-$$
-\min_x \left\Vert A x - b \right\Vert_2
-$$
+$QR = AP$ where $Q$ is an orthogonal matrix, $R$ is upper triangular, and $P$
+is a permutation matrix.
 @end tex
 @ifnottex
 
 @example
-min norm(A x - b)
+@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
-for overdetermined systems of equations (i.e.,
-@tex
-$A$
-@end tex
-@ifnottex
-@var{A}
-@end ifnottex
-is a tall, thin matrix).
 
-If only 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 are requested.)
+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 the matrix @var{A} is full, the permuted QR@tie{}factorization
-@code{[@var{Q}, @var{R}, @var{P}] = qr (@var{A})} forms the
-QR@tie{}factorization such that the diagonal entries of @var{R} are
-decreasing in magnitude order.  For example, given the matrix
-@code{a = [1, 2; 3, 4]},
+For example, given the matrix @code{@var{A} = [1, 2; 3, 4]},
 
 @example
 [@var{Q}, @var{R}, @var{P}] = qr (@var{A})
@@ -184,38 +197,68 @@
 @end group
 @end example
 
-The permuted @code{qr} factorization
-@code{[@var{Q}, @var{R}, @var{P}] = qr (@var{A})} factorization allows the
-construction of an orthogonal basis of @code{span (A)}.
-
-If the matrix @var{A} is sparse, then the sparse QR@tie{}factorization
-of @var{A} is computed using @sc{CSparse}.
-As the matrix @var{Q} is in general a full matrix,
-it is recommended to request only one return value,
-which is the @var{Q}-less factorization @var{R} of @var{A}, such that
-@code{@var{R} = chol (@var{A}' * @var{A})}.
+If the input matrix @var{A} is sparse then the sparse QR@tie{}factorization
+is computed using @sc{CSparse}.  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
+@var{R} such that @code{@var{R} = chol (@var{A}' * @var{A})}.
 
 If 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
+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})
-x = @var{R} \ @var{C}
+@var{x} = @var{R} \ @var{C}
 @end group
 @end example
 
-If the final argument is the scalar 0 and the number of rows is larger
-than the number of columns, then an @qcode{"economy"} factorization is
-returned, omitting zeroes of @var{R} and the corresponding columns of @var{Q}.
-That is, @var{R} will have only @code{size (@var{A},1)} rows.  In this case,
-@var{P} is a vector rather than a matrix.
+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 "economy" factorization is returned.
+When the original matrix @var{A} has size MxN and M > N then the "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 "economy"
+factorization the output @var{P} is always a vector rather than a matrix.
 
-If the final argument is the string @qcode{"vector"} then @var{P} is a
-permutation vector instead of a 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 */)
@@ -252,11 +295,12 @@
         }
       else if (args(nargin-1).is_string ())
         {
-          if (args(nargin-1).string_value () == "vector")
+          std::string str = args(nargin-1).string_value ();
+          if (str == "vector")
             vector_p = true;
-          else if (args(nargin-1).string_value () != "matrix")
+          else if (str != "matrix")
             error ("qr: type for P must be 'matrix' or 'vector', not %s",
-                   args(nargin-1).string_value ().c_str ());
+                   str.c_str ());
           have_b = (nargin > 2);
         }
       else if (! args(nargin-1).is_matrix_type ())
@@ -314,7 +358,7 @@
           if (arg.isreal ())
             {
               octave::math::qr<FloatMatrix>::type type
-                = qr_type<FloatMatrix> (nargin, nargout);
+                = qr_type<FloatMatrix> (nargout, economy);
 
               FloatMatrix m = arg.float_matrix_value ();
 
@@ -348,7 +392,7 @@
                   {
                     octave::math::qrp<FloatMatrix> fact (m, type);
 
-                    if (type == octave::math::qr<FloatMatrix>::economy || vector_p)
+                    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 ());
@@ -359,7 +403,7 @@
           else if (arg.iscomplex ())
             {
               octave::math::qr<FloatComplexMatrix>::type type
-                = qr_type<FloatComplexMatrix> (nargin, nargout);
+                = qr_type<FloatComplexMatrix> (nargout, economy);
 
               FloatComplexMatrix m = arg.float_complex_matrix_value ();
 
@@ -386,7 +430,7 @@
                 default:
                   {
                     octave::math::qrp<FloatComplexMatrix> fact (m, type);
-                    if (type == octave::math::qr<FloatComplexMatrix>::economy || vector_p)
+                    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 ());
@@ -399,7 +443,8 @@
         {
           if (arg.isreal ())
             {
-              octave::math::qr<Matrix>::type type = qr_type<Matrix> (nargin, nargout);
+              octave::math::qr<Matrix>::type type
+                = qr_type<Matrix> (nargout, economy);
 
               Matrix m = arg.matrix_value ();
 
@@ -432,7 +477,7 @@
                 default:
                   {
                     octave::math::qrp<Matrix> fact (m, type);
-                    if (type == octave::math::qr<Matrix>::economy || vector_p)
+                    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 ());
@@ -443,7 +488,7 @@
           else if (arg.iscomplex ())
             {
               octave::math::qr<ComplexMatrix>::type type
-                = qr_type<ComplexMatrix> (nargin, nargout);
+                = qr_type<ComplexMatrix> (nargout, economy);
 
               ComplexMatrix m = arg.complex_matrix_value ();
 
@@ -470,7 +515,7 @@
                 default:
                   {
                     octave::math::qrp<ComplexMatrix> fact (m, type);
-                    if (type == octave::math::qr<ComplexMatrix>::economy || vector_p)
+                    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 ());