dimforge
Condition number and nullspace functions
I am writing a piece of circuit simulation software that needs to detect ill-conditioned systems and calculate the nullspace of such a system (corresponding to invalid circuit nodes in my use case). Many linear algebra packages offer a nullspace function, such as Matlab's null(A, epsilon) and SciPy's null_space(A, relative_condition) which use the SVD to compute this, as well as a convenience method taking the ratio of the largest and smallest singular values to compute the condition number. I propose a pair of methods on the SVD<T, R, C> struct:
nullspace(&self, relative_condition) -> Matrix<T, R, D> where D is a dynamic dimension and the columnspace of the returned matrix corresponds to the nullspace of the original system.
condition_number(&self) -> T to get the condition number.
Adding this would also solve the root problem in #1125, since reduced row-echelon form is used to compute the nullspace in introductory linear algebra classes, but is numerically unstable and inconvenient for this on computers. This would also be a step toward the improved API discussed in #1117 .
Further API design considerations:
- Matlab uses an absolute epsilon, whereas SciPy uses a relative condition number. The second approach is probably better, but a variant that supports the same arguments as Matlab might be useful.
- Currently, the SVD struct optionally stores the right- and left- matricies in the decomposition, using Options. The suggested nullspace function would actually have to return an Option or Result; other methods use both, with no clear guidance on which is preferred. In the future, it would be nice to make this compile-time checked, since you probably know if you'll need the vectors when you compute the decomposition (and it's much cheaper to store them than compute them), but that's beyond the scope of this issue.
- Right now, LU decompositions are infallible, but may use a near-zero pivot, resulting in dubious output for ill-conditioned matricies. A version of the LU decomposition methods that reports their smallest pivot, or that take an epsilon and return Result, might be a logical counterpart to these methods on SVD.
I think both your suggestions for nullspace (or null_space? Bikeshedding time?) and condition_number are good, and I agree with your assessment that a relative criterion is more useful. If someone really wants an absolute value criterion then they can anyway implement this in terms of the relative version, as long as the docs are specific enough exactly what is being done.
However, some thought needs to be put into the three different cases of $m < n$, $m = n$, $m > n$ for a matrix $A \in \mathbb{R}^{m \times n}$, as this may possibly impact the design. I neither have the time nor the state of mind to think about this myself right now.
I agree that the SVD API is currently quite awkward due to the Option types used to indicate whether vectors are computed. On the other hand, doing this generically would complicate the SVD type quite a bit. A nice design here might nonetheless be worthwhile, but I suspect the "ideal" design might possibly include enums used in conjunction with const generics - but this is not yet a stable feature. It might be worth holding on to that before redesigning it - I am not sure.
As for LU decomposition being infallible - pivot value have no intuitive values known to me, at least, so it's very difficult to set anything reasonable here. At least it's completely unclear to me how pivot values are related to accuracy! Perhaps you have encountered something insightful on the topic? I think if you suspect that your linear system may be singular you rather ought to use something else than LU in the first place - such as SVD.
The good news about matrix sizing is that condition number calculation is independent of size (at least, so long as $\min(m, n) \geq 2$), since only the largest and smallest singular values are considered (and, notionally, represents how close a matrix or system is to having a rank less than the maximum permitted by dimension). The nullspace also should have a straightforward implementation, since it's just the vectors of the left ( $U$ ) matrix that have a near-zero corresponding entry on the diagonal of $\Sigma$, and you simply define that the diagonal entries with index above $\min(m, n)$ are zero.
As for the SVD API -- yes, or a struct SingularVectors and a struct SingularValuesOnly both impl SVDData, then adding that as a type parameter. Either way, a breaking and potentially messy change.
For LU decomposition -- The pivot values wind up on the diagonal of one of the output matrices (by convention, the U matrix). If a pivot is zero, then it will produce a zero on the diagonal of a triangular matrix, and that matrix cannot have full rank. So, the multiplied system cannot have full rank. This is not a failure per se of the LU decomposition, but it is useful information (a zero pivot is a sufficient, but not necessary, condition for a singular matrix) that might influence how much a user trusts the result. As you say, the LU decomposition shouldn't really be used for systems we expect to have a non-trivial nullspace, but it might be a useful diagnostic to expose.