Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

linalg: QR factorization #832

Merged
merged 35 commits into from
Sep 18, 2024
Merged

linalg: QR factorization #832

merged 35 commits into from
Sep 18, 2024

Conversation

perazz
Copy link
Contributor

@perazz perazz commented Jun 3, 2024

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 $.

  • base implementation
  • tests
  • exclude unsupported xdp
  • documentation
  • submodule
  • examples
  • all pure subroutine interfaces
  • 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
Copy link
Contributor Author

perazz commented Jun 4, 2024

From Fortran Monthly call:

  • 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 perazz marked this pull request as ready for review June 5, 2024 09:55
Copy link
Member

@jvdp1 jvdp1 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here are a few comments. LGTM. Thank you

doc/specs/stdlib_linalg.md Outdated Show resolved Hide resolved
doc/specs/stdlib_linalg.md Outdated Show resolved Hide resolved
doc/specs/stdlib_linalg.md Outdated Show resolved Hide resolved
doc/specs/stdlib_linalg.md Outdated Show resolved Hide resolved
doc/specs/stdlib_linalg.md Outdated Show resolved Hide resolved
doc/specs/stdlib_linalg.md Outdated Show resolved Hide resolved
doc/specs/stdlib_linalg.md Outdated Show resolved Hide resolved
doc/specs/stdlib_linalg.md Outdated Show resolved Hide resolved
src/stdlib_linalg.fypp Outdated Show resolved Hide resolved
src/stdlib_linalg_qr.fypp Outdated Show resolved Hide resolved
perazz and others added 7 commits July 8, 2024 08:55
Co-authored-by: Jeremie Vandenplas <[email protected]>
Co-authored-by: Jeremie Vandenplas <[email protected]>
Co-authored-by: Jeremie Vandenplas <[email protected]>
Co-authored-by: Jeremie Vandenplas <[email protected]>
Co-authored-by: Jeremie Vandenplas <[email protected]>
Co-authored-by: Jeremie Vandenplas <[email protected]>
Co-authored-by: Jeremie Vandenplas <[email protected]>
src/stdlib_linalg_qr.fypp Outdated Show resolved Hide resolved
test/linalg/test_linalg_qr.fypp Outdated Show resolved Hide resolved
@perazz
Copy link
Contributor Author

perazz commented Jul 11, 2024

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

Copy link
Member

@jvdp1 jvdp1 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

LGTM. Thank you @perazz for this addition.

@perazz
Copy link
Contributor Author

perazz commented Sep 15, 2024

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

@perazz
Copy link
Contributor Author

perazz commented Sep 18, 2024

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

@perazz perazz merged commit ccdba91 into fortran-lang:master Sep 18, 2024
16 checks passed
@perazz perazz deleted the qr branch September 18, 2024 06:20
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants