linalg: QR factorization
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 subroutineinterfaces - [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]):puresubroutine interfacecall qr_space(A, lwork [, err])query internal storage size for pre-allocation.
-
Special care was devoted to re-using internal storage such that a
pureand totally allocation-less method is available, if the user provides pre-allocated working arraystorage. Because LAPACK internals require two steps ("raw" factorization + ordered matrix output), temporary storage is borrowed from eithera,q, orrduring the operations. -
Both NumPy and SciPy have a cryptic UI that requires to provide a "method":
full,reducedraweconomic,r. I believe this is hard to understand. Instead, the current version proposes to decide it autonomously:
- If
shape(Q)==[m,m]andshape(R)==[m,n], perform a full factorization, $A=QR$. - If
shape(Q)==[m,k]andshape(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
From Fortran Monthly call:
- [x] require that
QandRhave exact sizes for either the full or the reduced problem (do not just check that there is "enough" storage)
Thank you for your edits @jvdp1 @jalvesz, looks much better now.
Thank you @jvdp1. I think it makes sense to wait for a bit longer, and then merge if there are no further comments.
Thank you @jvdp1 @jalvesz, I will merge due to no outstanding comments.