Suppose that we don't know $n \times n$ matrix $A$ explicitly but we are only able to compute products $Ax$ where $x$ is a column vector with $n$ elements. Is there an algorithm to determine whether $A$ is singular?
-
4$\begingroup$ There is, of course, the obvious algorithm that computes the products of $A$ with the $n$ columns of the identity matrix (canonical basis) and forms the matrix explicitly. Do you have constraints that disallow it? $\endgroup$– Federico PoloniCommented Aug 4, 2020 at 11:56
-
$\begingroup$ I'd like to suggest working the requirement about only knowing products Ax into the title somehow. As it stands, the title suggests a more generic problem. $\endgroup$– Robert DodierCommented Aug 5, 2020 at 5:03
-
$\begingroup$ When you have a concrete problem, you can look if there are "rank revealing" algorithms for solving it. $\endgroup$– alloCommented Aug 5, 2020 at 10:17
-
$\begingroup$ Not enough reputation to make a comment. Can you identify column A_i by setting x = e_i (basis vector i)? Once you have the whole matrix A just evaluate it's eigenvalues. $\endgroup$– BalavsCommented Aug 5, 2020 at 17:07
4 Answers
If you can compute products with $A$ and $A^T$, as you specify in a comment, you can run the classical sparse SVD algorithms such as scipy.sparse.linalg.svds
, Matlab's svds
, or Julia's Arpack.svds
, which are based on Lanczos bidiagonalization. They are designed to compute singular values, and are likely to be more robust than a minimization routine coded by hand.
Then the distance (in Euclidean or Frobenius norm) from $A$ to the nearest singular matrix is precisely the smallest singular value $\sigma_{\min}$. You won't be able to return an exact yes/no answer using IEEE arithmetic, so this is the best you can do.
-
$\begingroup$ Interval or radial arithmetic with outward rounding ,might be able to definitively determine a matrix iis not singular. But I don't believe it would be able to definitivlely determine it is singular. I.e., two possible results: 1) not singular 2) may or may not be sngular. See for example ti3.tu-harburg.de/paper/rump/Ru10a_alt.pdf using INTLAB under MATLAB. for verified computation of singular values. $\endgroup$ Commented Aug 4, 2020 at 17:52
-
$\begingroup$ @MarkL.Stone I consider "interval arithmetic" to be something distinct from "IEEE arithmetic", even if it is implemented using it. $\endgroup$ Commented Aug 4, 2020 at 19:11
-
$\begingroup$ Yes, but I was giving an alternative to IEEE arithmetic. The OP did not state that only IEEE arithmetic was allowed. $\endgroup$ Commented Aug 5, 2020 at 2:50
-
$\begingroup$ I have now provided an expanded version of my comment as an answer. $\endgroup$ Commented Aug 5, 2020 at 10:52
-
$\begingroup$ I solved the problem by computing the eigenvalue with smallest magnitude with ARPACK. Is this correct? $\endgroup$– tohoynCommented Aug 5, 2020 at 11:09
I also suggest looking into the condition number estimators, which will (with some degree of [un]reliability) predict how effectively numerically singular the matrix is.
In particular,
- "Spectral Condition-Number Estimation of Large Sparse Matrices" based on LSQR seems like an interesting choice.
attracted my attention one day. I would also suggest going over the references in this paper to get some perspective on the problem.
@Federico Poloni 's fine answer states the impossibility of getting an exact yes/no answer using IEEE arithmetic.
However, using interval arithmetic with outward rounding, it is possible to get a "not singular/don't know" answer. In particular, it may be possible to definitively conclude that the smallest singular value is strictly greater than zero.
This interval arithmetic computation using outward rounding of the smallest singular value, can be done, for example, using INTLAB (developed by Siegfried M. Rump) under MATLAB.
The result will be either:
- "not singular", i.e., the matrix is not singular, because its smallest singular value is definitely > 0
or
- "don't know", i.e., not definitively determined whether or not the matrix is singular, because not definitively determined whether or not its smallest singular value = 0
"Standard" interval arithmetic with outward rounding (to include, INTLAB) is done using double precision floating point arithmetic. However, it is possible to perform interval arithmetic with outward rounding in higher precision. Doing so might allow a definitive "not singular" result for some matrices which "evaluate" to "don't know" using double precision interval arithmetic with outward rounding. In the (unimplementable) limit of infinite precision interval arithmetic with outward rounding, it can be definitively determined whether or not a given matrix is singular.
Here's my 2 cents. I would set up the following minimization problem
$$ \pi(x) = \frac{1}{2} (Ax)^T(Ax) $$
If $A$ has eigenvalues which are zero, there will exist a nonzero $x$ such that $\pi(x)=0$. So, I would try computing the gradient of $\pi$ wrt $x$ and use a gradient descent algorithm to drive $\pi$ towards zero. If you get reasonably close to zero ($\pi\approx$ 1e-12), then the matrix is singular. The first variation of $\pi$ can be computed to be
$$ \delta\pi = x^TA^TA\delta{x} = (Ax)^TA\delta{x} = g^T\delta{x}, $$ where $g$ is the gradient. So $g$ is $$ g = A^TAx $$
You'd also need to avoid the $x=0$ case. Starting from a non zero random vector might help.
Here, you need to compute the product of $A^T$ with $Ax$. Given the constraint that you have ( you can only compute $Ax$), I'm not sure if this is possible. Maybe the experts can answer.
-
$\begingroup$ Thanks for your answer. I can actually compute $A^T x$, too. $\endgroup$– tohoynCommented Aug 4, 2020 at 10:48
-
$\begingroup$ Won’t this only work if you maintain a fixed (positive) norm for $x$? (Otherwise $x=0$ always works.) This constraint makes the optimization problem quite a bit harder. $\endgroup$– cdipaoloCommented Aug 4, 2020 at 17:57
-
$\begingroup$ @cdipaolo You're absolutely right. Maybe an optimization algorithm capable of handling the constraint $x^Tx=1$. But Frederico-Polini's answer is the best solution. $\endgroup$– NNNCommented Aug 5, 2020 at 6:31
-
$\begingroup$ If you determine a singular value $< 10^{-12}$ in $A^TA$, then that only gives a singular value $< 10^{-6}$ in $A$. Depending on the situation, this may be too large to conclude singularity. $\endgroup$ Commented Aug 5, 2020 at 7:40
-
$\begingroup$ You might also want to try SlepC: slepc.upv.es/documentation/current/docs/manualpages/EPS/… $\endgroup$– NNNCommented Aug 5, 2020 at 12:00