From bug-octave-request at bevo dot che dot wisc dot edu Fri Dec 14 12:28:55 2001
Subject: economy QR isn't economical
From: Paul Kienzle <pkienzle at jazz dot ncnr dot nist dot gov>
To: bug-octave at bevo dot che dot wisc dot edu
Date: Fri, 14 Dec 2001 13:28:37 -0500

Please find enclosed a patch to the liboctave QR routines so that they only
generate the requested columns of the QR decomposition.  In addition, the
*QRP functions are revised to more closely match the structure of the *QR
functions.

I have tested the patch on real and complex inputs with m>n and m<n with
and without permutations for both the economy and the full cases.  I did
not test the "raw" matrix returned by R=qr(A).

function testqr(q,r,a,p)
  if (rows(a)<columns(a)) s='wide '; else s='tall '; endif
  if (isreal(a)) t='real '; else t='complex '; endif
  if (nargin == 3)
     disp(['straight ', s, t]);
     printf("|q*r-a|      =%g\n",norm(q*r-a));
     printf("|q'*q-I|     =%g\n",norm(q'*q-eye(columns(q))));
  else
     disp(['permuted ', s, t]);
     if (is_vector(p))
        printf("|q*r-a(:,p)| =%g\n",norm(q*r-a(:,p)));
     else
	printf("|q*r-a*p|    =%g\n", norm(q*r - a*p));
     endif
     printf("|q'*q-I|     =%g\n", norm(q'*q-eye(columns(q))));
  endif
  fflush(stdout);
endfunction

disp('economy');
a = rand(5000,20);
[q,r]=qr(a,0); testqr(q,r,a);
[q,r]=qr(a',0); testqr(q,r,a');
[q,r,p]=qr(a,0); testqr(q,r,a,p);
[q,r,p]=qr(a',0); testqr(q,r,a',p);

a = a+1i*randn(size(a))*sqrt(eps);
[q,r]=qr(a,0); testqr(q,r,a);
[q,r]=qr(a',0); testqr(q,r,a');
[q,r,p]=qr(a,0); testqr(q,r,a,p);
[q,r,p]=qr(a',0); testqr(q,r,a',p);

disp('full');
a = [ ones(1,15); sqrt(eps)*eye(15) ];
[q,r]=qr(a); testqr(q,r,a);
[q,r]=qr(a'); testqr(q,r,a');
[q,r,p]=qr(a); testqr(q,r,a,p);
[q,r,p]=qr(a'); testqr(q,r,a',p);

a = a+1i*randn(size(a))*sqrt(eps);
[q,r]=qr(a); testqr(q,r,a);
[q,r]=qr(a'); testqr(q,r,a');
[q,r,p]=qr(a); testqr(q,r,a,p);
[q,r,p]=qr(a'); testqr(q,r,a',p);

Paul Kienzle
pkienzle at users dot sf dot net


Index: CmplxQR.cc
===================================================================
RCS file: /cvs/octave/liboctave/CmplxQR.cc,v
retrieving revision 1.21
diff -c -p -r1.21 CmplxQR.cc
*** CmplxQR.cc	1997/07/10 23:34:06	1.21
--- CmplxQR.cc	2001/12/14 18:13:00
*************** ComplexQR::init (const ComplexMatrix& a,
*** 74,80 ****
    int info = 0;
  
    ComplexMatrix A_fact;
!   if (m > n)
      {
        A_fact.resize (m, m);
        A_fact.insert (a, 0, 0);
--- 74,80 ----
    int info = 0;
  
    ComplexMatrix A_fact;
!   if (m > n && qr_type != QR::economy)
      {
        A_fact.resize (m, m);
        A_fact.insert (a, 0, 0);
*************** ComplexQR::init (const ComplexMatrix& a,
*** 106,123 ****
  	}
        else
  	{
! 	  volatile int n2;
  
  	  if (qr_type == QR::economy && m > n)
! 	    {
! 	      n2 = n;
! 	      r.resize (n, n, 0.0);
! 	    }
  	  else
! 	    {
! 	      n2 = m;
! 	      r.resize (m, n, 0.0);
! 	    }
  
  	  for (int j = 0; j < n; j++)
  	    {
--- 106,120 ----
  	}
        else
  	{
! 	  int n2 = m;
  
+ 	  if (qr_type == QR::economy)
+ 	    n2 = min_mn;
+ 
  	  if (qr_type == QR::economy && m > n)
! 	    r.resize (n, n, 0.0);
  	  else
! 	    r.resize (m, n, 0.0);
  
  	  for (int j = 0; j < n; j++)
  	    {
*************** ComplexQR::init (const ComplexMatrix& a,
*** 126,136 ****
  		r.elem (i, j) = A_fact.elem (i, j);
  	    }
  
! 	  lwork = 32*m;
  	  work.resize (lwork);
  	  Complex *pwork = work.fortran_vec ();
  
! 	  F77_XFCN (zungqr, ZUNGQR, (m, m, min_mn, tmp_data, m, ptau,
  				     pwork, lwork, info));
  
  	  if (f77_exception_encountered)
--- 123,133 ----
  		r.elem (i, j) = A_fact.elem (i, j);
  	    }
  
! 	  lwork = 32*n2;
  	  work.resize (lwork);
  	  Complex *pwork = work.fortran_vec ();
  
! 	  F77_XFCN (zungqr, ZUNGQR, (m, n2, min_mn, tmp_data, m, ptau,
  				     pwork, lwork, info));
  
  	  if (f77_exception_encountered)
Index: CmplxQRP.cc
===================================================================
RCS file: /cvs/octave/liboctave/CmplxQRP.cc,v
retrieving revision 1.20
diff -c -p -r1.20 CmplxQRP.cc
*** CmplxQRP.cc	1997/03/27 16:18:34	1.20
--- CmplxQRP.cc	2001/12/14 18:13:00
*************** ComplexQRP::init (const ComplexMatrix& a
*** 68,74 ****
        return;
      }
  
!   Array<Complex> tau (m < n ? m : n);
    Complex *ptau = tau.fortran_vec ();
  
    int lwork = 3*n > 32*m ? 3*n : 32*m;
--- 68,75 ----
        return;
      }
  
!   int min_mn = m < n ? m : n;
!   Array<Complex> tau (min_mn);
    Complex *ptau = tau.fortran_vec ();
  
    int lwork = 3*n > 32*m ? 3*n : 32*m;
*************** ComplexQRP::init (const ComplexMatrix& a
*** 78,84 ****
    int info = 0;
  
    ComplexMatrix A_fact = a;
!   if (m > n)
      A_fact.resize (m, m, 0.0);
  
    Complex *tmp_data = A_fact.fortran_vec ();
--- 79,85 ----
    int info = 0;
  
    ComplexMatrix A_fact = a;
!   if (m > n && qr_type != QR::economy)
      A_fact.resize (m, m, 0.0);
  
    Complex *tmp_data = A_fact.fortran_vec ();
*************** ComplexQRP::init (const ComplexMatrix& a
*** 114,136 ****
  	    p.elem (jpvt.elem (j) - 1, j) = 1.0;
  	}
  
        if (qr_type == QR::economy && m > n)
  	r.resize (n, n, 0.0);
        else
  	r.resize (m, n, 0.0);
  
-       int min_mn = m < n ? m : n;
- 
        for (int j = 0; j < n; j++)
  	{
  	  int limit = j < min_mn-1 ? j : min_mn-1;
  	  for (int i = 0; i <= limit; i++)
  	    r.elem (i, j) = A_fact.elem (i, j);
  	}
- 
-       int n2 = m;
-       if (qr_type == QR::economy)
- 	n2 = min_mn;
  
        F77_XFCN (zungqr, ZUNGQR, (m, n2, min_mn, tmp_data, m, ptau,
  				 pwork, lwork, info));
--- 115,135 ----
  	    p.elem (jpvt.elem (j) - 1, j) = 1.0;
  	}
  
+       int n2 = m;
+       if (qr_type == QR::economy)
+ 	n2 = min_mn;
+ 
        if (qr_type == QR::economy && m > n)
  	r.resize (n, n, 0.0);
        else
  	r.resize (m, n, 0.0);
  
        for (int j = 0; j < n; j++)
  	{
  	  int limit = j < min_mn-1 ? j : min_mn-1;
  	  for (int i = 0; i <= limit; i++)
  	    r.elem (i, j) = A_fact.elem (i, j);
  	}
  
        F77_XFCN (zungqr, ZUNGQR, (m, n2, min_mn, tmp_data, m, ptau,
  				 pwork, lwork, info));
Index: dbleQR.cc
===================================================================
RCS file: /cvs/octave/liboctave/dbleQR.cc,v
retrieving revision 1.21
diff -c -p -r1.21 dbleQR.cc
*** dbleQR.cc	1997/07/10 23:34:06	1.21
--- dbleQR.cc	2001/12/14 18:13:00
*************** QR::init (const Matrix& a, QR::type qr_t
*** 71,84 ****
  
    int info = 0;
  
!   Matrix A_fact;
!   if (m > n)
!     {
!       A_fact.resize (m, m);
!       A_fact.insert (a, 0, 0);
!     }
!   else
!     A_fact = a;
  
    double *tmp_data = A_fact.fortran_vec ();
  
--- 71,79 ----
  
    int info = 0;
  
!   Matrix A_fact = a;
!   if (m > n && qr_type != QR::economy)
!       A_fact.resize (m, m, 0.0);
  
    double *tmp_data = A_fact.fortran_vec ();
  
*************** QR::init (const Matrix& a, QR::type qr_t
*** 104,134 ****
  	}
        else
  	{
! 	  volatile int n2;
  
  	  if (qr_type == QR::economy && m > n)
! 	    {
! 	      n2 = n;
! 	      r.resize (n, n, 0.0);
! 	    }
  	  else
! 	    {
! 	      n2 = m;
! 	      r.resize (m, n, 0.0);
! 	    }
  
  	  for (int j = 0; j < n; j++)
  	    {
  	      int limit = j < min_mn-1 ? j : min_mn-1;
  	      for (int i = 0; i <= limit; i++)
! 		r.elem (i, j) = tmp_data[m*j+i];
  	    }
  
! 	  lwork = 32*m;
  	  work.resize (lwork);
  	  double *pwork = work.fortran_vec ();
  
! 	  F77_XFCN (dorgqr, DORGQR, (m, m, min_mn, tmp_data, m, ptau,
  				     pwork, lwork, info));
  
  	  if (f77_exception_encountered)
--- 99,125 ----
  	}
        else
  	{
! 	  int n2 = m;
! 	  if (qr_type == QR::economy)
! 	    n2 = min_mn;
  
  	  if (qr_type == QR::economy && m > n)
! 	    r.resize (n, n, 0.0);
  	  else
! 	    r.resize (m, n, 0.0);
  
  	  for (int j = 0; j < n; j++)
  	    {
  	      int limit = j < min_mn-1 ? j : min_mn-1;
  	      for (int i = 0; i <= limit; i++)
! 		r.elem (i, j) = A_fact.elem(i, j);
  	    }
  
! 	  lwork = 32*n2;
  	  work.resize (lwork);
  	  double *pwork = work.fortran_vec ();
  
! 	  F77_XFCN (dorgqr, DORGQR, (m, n2, min_mn, tmp_data, m, ptau,
  				     pwork, lwork, info));
  
  	  if (f77_exception_encountered)
Index: dbleQRP.cc
===================================================================
RCS file: /cvs/octave/liboctave/dbleQRP.cc,v
retrieving revision 1.20
diff -c -p -r1.20 dbleQRP.cc
*** dbleQRP.cc	1997/03/27 16:18:46	1.20
--- dbleQRP.cc	2001/12/14 18:13:00
*************** QRP::init (const Matrix& a, QR::type qr_
*** 67,73 ****
        return;
      }
  
!   Array<double> tau (m < n ? m : n);
    double *ptau = tau.fortran_vec ();
  
    int lwork = 3*n > 32*m ? 3*n : 32*m;
--- 67,74 ----
        return;
      }
  
!   int min_mn = m < n ? m : n;
!   Array<double> tau (min_mn);
    double *ptau = tau.fortran_vec ();
  
    int lwork = 3*n > 32*m ? 3*n : 32*m;
*************** QRP::init (const Matrix& a, QR::type qr_
*** 77,83 ****
    int info = 0;
  
    Matrix A_fact = a;
!   if (m > n)
      A_fact.resize (m, m, 0.0);
  
    double *tmp_data = A_fact.fortran_vec ();
--- 78,84 ----
    int info = 0;
  
    Matrix A_fact = a;
!   if (m > n && qr_type != QR::economy)
      A_fact.resize (m, m, 0.0);
  
    double *tmp_data = A_fact.fortran_vec ();
*************** QRP::init (const Matrix& a, QR::type qr_
*** 109,131 ****
  	    p.elem (jpvt.elem (j) - 1, j) = 1.0;
  	}
  
        if (qr_type == QR::economy && m > n)
  	r.resize (n, n, 0.0);
        else
  	r.resize (m, n, 0.0);
  
-       int min_mn = m < n ? m : n;
- 
        for (int j = 0; j < n; j++)
  	{
  	  int limit = j < min_mn-1 ? j : min_mn-1;
  	  for (int i = 0; i <= limit; i++)
  	    r.elem (i, j) = A_fact.elem (i, j);
  	}
- 
-       int n2 = m;
-       if (qr_type == QR::economy)
- 	n2 = min_mn;
  
        F77_XFCN (dorgqr, DORGQR, (m, n2, min_mn, tmp_data, m, ptau,
  				 pwork, lwork, info));
--- 110,130 ----
  	    p.elem (jpvt.elem (j) - 1, j) = 1.0;
  	}
  
+       int n2 = m;
+       if (qr_type == QR::economy)
+ 	n2 = min_mn;
+ 
        if (qr_type == QR::economy && m > n)
  	r.resize (n, n, 0.0);
        else
  	r.resize (m, n, 0.0);
  
        for (int j = 0; j < n; j++)
  	{
  	  int limit = j < min_mn-1 ? j : min_mn-1;
  	  for (int i = 0; i <= limit; i++)
  	    r.elem (i, j) = A_fact.elem (i, j);
  	}
  
        F77_XFCN (dorgqr, DORGQR, (m, n2, min_mn, tmp_data, m, ptau,
  				 pwork, lwork, info));



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------