stdlib icon indicating copy to clipboard operation
stdlib copied to clipboard

linalg: QR factorization

Open perazz opened this issue 1 year ago • 2 comments

Compute the QR factorization of a real or complex matrix: $A = Q R $, where q is orthonormal and r is upper-triangular. Matrix $A$ has size [m,n], with $m\ge n$. Based on LAPACK General QR factorization (*GEQRF) and ordered matrix output (*ORGQR, *UNGQR). Option for full or reduced factorization: given k = min(m,n), one can write $A = ( Q_1 Q_2 ) \cdot ( \frac{R_1}{0})$. The user may want the full problem or the reduced problem only $A = Q_1 R_1 $.

  • [x] base implementation
  • [x] tests
  • [x] ~~exclude unsupported xdp~~
  • [x] documentation
  • [x] submodule
  • [x] examples
  • [x] all pure subroutine interfaces
  • [x] workspace size evaluation

Prior art

  • Numpy: linalg.qr(a, mode={'reduced', 'complete', 'r', 'raw')
  • Scipy: scipy.linalg.qr(a, overwrite_a=False, lwork=None, mode='full', pivoting=False, check_finite=True)

Proposed implementation

  • call qr(A,Q,R [, overwrite_a] [, storage] [, err]) : pure subroutine interface
  • call qr_space(A, lwork [, err]) query internal storage size for pre-allocation.
  1. Special care was devoted to re-using internal storage such that a pure and totally allocation-less method is available, if the user provides pre-allocated working array storage. Because LAPACK internals require two steps ("raw" factorization + ordered matrix output), temporary storage is borrowed from either a, q, or r during the operations.

  2. Both NumPy and SciPy have a cryptic UI that requires to provide a "method": full, reduced raw economic, r. I believe this is hard to understand. Instead, the current version proposes to decide it autonomously:

  • If shape(Q)==[m,m] and shape(R)==[m,n], perform a full factorization, $A=QR$.
  • If shape(Q)==[m,k] and shape(R)==[k,n], perform a "reduced" factorization, $A=Q_1R_1$.
  • On insufficient/unfit size(s), raise an error.

In other words, the user decides what method is required based on the size of the arrays passed to qr.

cc: @fortran-lang/stdlib @jvdp1 @jalvesz @everythingfunctional

perazz avatar Jun 03 '24 08:06 perazz

From Fortran Monthly call:

  • [x] require that Q and R have exact sizes for either the full or the reduced problem (do not just check that there is "enough" storage)

perazz avatar Jun 04 '24 17:06 perazz

Thank you for your edits @jvdp1 @jalvesz, looks much better now.

perazz avatar Jul 11 '24 13:07 perazz

Thank you @jvdp1. I think it makes sense to wait for a bit longer, and then merge if there are no further comments.

perazz avatar Sep 15 '24 12:09 perazz

Thank you @jvdp1 @jalvesz, I will merge due to no outstanding comments.

perazz avatar Sep 18 '24 06:09 perazz