From: hui@yrloc.ipsa.reuter.COM (Roger Hui)
Subject: Re: Numerical Analysis + APL/J
Organization: Iverson Software Inc.
Date: Fri, 14 Aug 92 03:26:19 GMT


Roger Hui writes:
   My msg said, "QR applies to a real or complex matrix A with linearly-
   independent columns and returns an Hermitian matrix Q and an upper
   triangular matrix R, individually boxed, such that A = Q +/ .* R."
   I would be quite surprised if this can be done in O(n^2) time.
 
Andrew Mullhaupt writes:
   It can't. When you claim that this can benefit from Strassen
   multiplication, I assumed that you were trying to speed up the
   application of the individual steps. (Which are examples of O(n^2)
   matrix multiplications which cannot be helped by Strassen...)
 
   The QR factorization is O(n^3), (actually O(mn^2)), and as far as
   I can tell there is _nothing_ in the algorithm which can benefit
   from Strassen's trick.
 
   What exactly are you claiming Strassen multiplication, or something
   like it, can do for the QR factorization?
 
[Aside:  You say, "It can't".  I believe it would be quite difficult to
prove that QR can't be done in time O(n^2) (or O(m*n) on an (m,n) matrix).]
 
In my original msg, I said (and you quoted), "The algorithms show
that RI, QR, LU, and Cholesky are of the same (sequential) complexity
as matrix multiplication and therefore can benefit directly from
Strassen and other fast matrix multiplications."  I meant exactly that:
If two n,n matrices can be multiplied in time O(n^e), then QR on an m,n
matrix takes time O(m*n^e-1).  If m=n, QR take time O(n*n^e-1)=O(n^e),
the same as matrix multiplication.  For example, if matrices are multiplied
using the Strassen algorithm, then QR takes time O(m*n^(log 7)-1).
(log=.2&^. , i.e. base-2 log.)
 
To substantiate my claim, I've reproduced here QR from my original msg,
annotated with the required times to within the same order.  (To avoid
confusion, I've renamed a local variable in q2 from m to n2.)  m T n is
the time required by QR on an m,n matrix.
 
ip =. +/ .*
 
t=.0 0$''
t=.t, 'n  =. en y.'                    NB.  1
t=.t, 'n2 =. >.-: n'                   NB.  1
t=.t, 'a0 =. n2{."1 y.'                NB.  m*n
t=.t, 'a1 =. n2}."1 y.'                NB.  m*n
t=.t, 't0 =. QR a0'                    NB.  m T n%2
t=.t, 'q0 =. >@{. t0'                  NB.  m*n
t=.t, 'r0 =. >@{: t0'                  NB.  m*n
t=.t, 'c  =. (+|:q0) ip a1'            NB.  (m*n)+ (m%n%2)*(n%2)^e
t=.t, 't1 =. QR a1 - q0 ip c'          NB.  (m*n)+((m%n%2)*(n%2)^e)+m T n%2
t=.t, 'q1 =. >@{. t1'                  NB.  m*n
t=.t, 'r1 =. >@{: t1'                  NB.  m*n
t=.t, '(q0,.q1);(r0,.c),(-n){."1 r1'   NB.  m*n
q2   =. t : ''
 
en   =. 1&{@(,&1 1)@$                  NB.  1
norm =. (%:@ip +)@,                    NB.  m
QR   =. q2`((% ;&,. ,~@en@[ $ ]) norm) @. (1&>:@en)
 
Therefore, m T n <-> (2*m T n%2) + (2*(m%n%2)*(n%2)^e) + O(m*n) and
m T 1 <-> O(m), and an inductive argument shows that T is O(m*n^e-1).
[I assume e>2; if e=2, other terms dominate, and T is O(m*n*log n).]
 
Andrew, despite the reservations you'd previously expressed in this forum,
I'd like to suggest that it may be worthwhile for you to get a copy of J,
if only to use as a tool of thought in verifying the workings of QR.
 
p.s. So how _do_ you compute the primes up to n in less than O(n) time?

------------------------------------
Roger Hui, Iverson Software Inc., 33 Major Street, Toronto, Ontario  M5S 2K9
Phone: (416) 925 6096;  Fax: (416) 488 7559