Skip to content

Some linalg functions are not well-defined #207

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

Closed
lezcano opened this issue Jun 29, 2021 · 2 comments · Fixed by #224
Closed

Some linalg functions are not well-defined #207

lezcano opened this issue Jun 29, 2021 · 2 comments · Fixed by #224
Labels
topic: Linear Algebra Linear algebra.
Milestone

Comments

@lezcano
Copy link
Contributor

lezcano commented Jun 29, 2021

The functions eigh, svd, qr and lstsq are not well-defined as written.

Uniqueness lstsq

linalg.lstsq docs currently read:
“Returns the least-squares solution to a linear matrix equation”
but the solution needn't be unique whenever A has non-trivial kernel (e.g. A = 0).

Proposed solution:

  • Rewrite it as "Returns a least-square solution [...]"
  • Perhaps we want to add a note on how to handle these cases (e.g. the solution with smallest norm is chosen) so that this function is well-defined for every input matrix. Note that this would not allow for the implementation of this function via faster but less robust algorithms such as LAPACK's gels, which assumes that A is full-rank and performs a QR decomposition.

Uniqueness eigh and svd:

  • Eigenvectors (resp. singular vectors) of a generic matrix are only defined up to sign in the real case and multiplication by a scalar of norm 1 in the complex case
  • When the matrix has repeated eigenvalues (resp. singular values) the associated eigenvectors are just defined up to multiplication by a rotation (or multiplication by a unitary matrix in the complex case).

Proposed solution:

  • Change the wording from "Returns the eigenvalues and eigenvectors of a symmetric matrix" to "Returns an eigendecomposition of a symmetric matrix".
  • Change the wording of linalg.svd as well.

QR decomposition cannot be defined

The qr decomposition is just well-defined whenever the matrix has full column rank. In the singular case pivoting is necessary.

Proposed solution: The behaviour of this function should be marked as unspecified whenever the input matrix does not have full column rank.

cc @rgommers @mruberry

@rgommers rgommers added the topic: Linear Algebra Linear algebra. label Jun 30, 2021
@rgommers
Copy link
Member

Thanks for the detailed review @lezcano. Your proposed solutions for eigh, svd and qr all sound good to me, as does the first bullet for lstsq.

Perhaps we want to add a note on how to handle these cases (e.g. the solution with smallest norm is chosen) so that this function is well-defined for every input matrix. Note that this would not allow for the implementation of this function via faster but less robust algorithms such as LAPACK's gels, which assumes that A is full-rank and performs a QR decomposition.

I don't think we should add such a note. One principle we'd tried to adhere to in the API design is "don't make decisions that limit implementers in how they optimize for performance". So taking a decision on performance-robustness tradeoffs for lstsq doesn't sound optimal.

Pointing out the issue is good though - but I'd do it as a note like "if the solution is not unique, implementations may return different results".

@lezcano
Copy link
Contributor Author

lezcano commented Jul 1, 2021

For what is worth, I do not know of any implementation of lstsq that does not return the minimum norm solution. For example, all the LAPACK implementations do so. Note that even gels falls in this category as well, albeit under the assumption that the matrix is full-rank.

On a similar decision to this one, see the OP of the PR for linalg.matrix_rank. The OP says that, even though some implementations provide an hermitian boolean flag, this is not included in the specification because it denotes an optional optimisation. In the same way, the behaviour of gels could also be provided by the libraries via an optional flag full_rank=False or similar.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
topic: Linear Algebra Linear algebra.
Projects
None yet
Development

Successfully merging a pull request may close this issue.

3 participants