8
$\begingroup$

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?

$\endgroup$
  • 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 Poloni Aug 4 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 Dodier Aug 5 at 5:03
  • $\begingroup$ When you have a concrete problem, you can look if there are "rank revealing" algorithms for solving it. $\endgroup$ – allo Aug 5 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$ – Balavs Aug 5 at 17:07
12
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
  • $\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$ – Mark L. Stone Aug 4 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$ – Federico Poloni Aug 4 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$ – Mark L. Stone Aug 5 at 2:50
  • $\begingroup$ I have now provided an expanded version of my comment as an answer. $\endgroup$ – Mark L. Stone Aug 5 at 10:52
  • $\begingroup$ I need to do the computations in C++ and the matrices are not given explicitly. I'm doing ARPACK iteration for vector $(x_n^1,...,x_n^N,y_n^1,...,y_n^N)$ with $$\mathbf{x}_{n+1} = A \mathbf{y}_n$$ and $$\mathbf{y}_{n+1} = A^T \mathbf{x}_n$$ I set the sigma parameter to "SM" in order to compute the smallest (absolute value) singular values. Is this correct? $\endgroup$ – Tommi Höynälänmaa Aug 6 at 5:38
2
$\begingroup$

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,

attracted my attention one day. I would also suggest going over the references in this paper to get some perspective on the problem.

| cite | improve this answer | |
$\endgroup$
2
$\begingroup$

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

See Verified bounds for singular values, in particular for the spectral norm of a matrix and its inverse, Siegfried M. Rump, BIT, 51(2):367384, 2011.

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:

  1. "not singular", i.e., the matrix is not singular, because its smallest singular value is definitely > 0

or

  1. "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.

| cite | improve this answer | |
$\endgroup$
1
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thanks for your answer. I can actually compute $A^T x$, too. $\endgroup$ – Tommi Höynälänmaa Aug 4 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$ – cdipaolo Aug 4 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$ – Nachiket Aug 5 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$ – Lutz Lehmann Aug 5 at 7:40
  • $\begingroup$ You might also want to try SlepC: slepc.upv.es/documentation/current/docs/manualpages/EPS/… $\endgroup$ – Nachiket Aug 5 at 12:00
0
$\begingroup$

I solved the problem by computing the eigenvalue with smallest magnitude with ARPACK. Is this correct?

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Sort of. The smallest eigenvalue does not provide a meaningful metric of how close a matrix is to singularity, unlike the smallest singular value. $\endgroup$ – Federico Poloni Aug 5 at 11:44
  • $\begingroup$ This does not provide an answer to the question. To critique or request clarification from an author, leave a comment below their post. - From Review $\endgroup$ – Alone Programmer Aug 7 at 16:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.