1

Dear numerical algebra experts,

I am trying to find $\alpha$ for which $A(\alpha)x=0$ has a non-trivial solution ($x\neq0$), i.e. I am looking for $\alpha_0$ for which the null space of $A(\alpha_0)$ is not empty. A (square, NxN, real) is the coefficient matrix of a physics problem (constructed from boundary conditions and an Ansatz that solves the corresponding diff. equations of interest). My general strategy (performed in Numpy/Scipy but would certainly work analogously in different environments) is

  • search for the root $\alpha_0$ in the determinant of $A(\alpha)$
  • Perform a SVD of $A(\alpha_0)$: $A(\alpha_0) = UsV$ and compile the null space from $V$ (the row vector that corresponds to the smallest singular value $\sigma$ in $s$).

This strategy works well for certain parameter spaces of $\alpha$: the final result is what I expect (physics-wise) and all numerical parameters behave well.

However, for other parameter spaces of $\alpha$, I observe that I can still find $\alpha_0$ (the root of the determinant) but the final result (the "null space" given by the last row vector of V) does not behave well anymore (discontinuities). After thoroughly checking the behavior of the relevant numerical figure of merits, I figured out that this problem seems to be related to $A(\alpha_0)$ becoming rank-deficient (in a numerical sense): using a common threshold definition of max($\sigma)\times N \times$ (machine epsilon), the described problem starts to occur once the second (or even third) smallest $\sigma$ exceeds this threshold by a certain value (cannot exactly nail it).

In conclusion, it seems there are issues with the floating-point precision during the SVD. Since working on a machine with better precision is not really on option, I am somehow stuck at this point. Is there any strategy (or hint) on how to proceed?

Any feedback is highly appreciated!

Update:

I followed the nice idea proposed by Ben Grossmann. However, it did not work out completely since the scaling of $A(\alpha)$ being ill-conditioned seems to be highly unfortunate. Nevertheless, I figured out that some elements of $A$ scaled (in magnitude) signifianctly worse than the others (which might explain why a global rescaling/preconditioning did not fully work). So I did a manual rescaling (preconditioning) of only these elements for an $\alpha$ in the center of the parameter space of interest. This worked out quite nicely and allowed me to solve the issue. Thanks @ Ben.

Michael Steinke
  • 11
  • 2
  • Please clarify what exactly $A(\alpha)$ is supposed to be. My best guess is that $A:\Bbb R \to \Bbb R^{N \times N}$ is a continuous function (i.e. a matrix whose entries depend on an input parameter), so that for $\alpha \in \Bbb R$, $A(\alpha)$ is an $N \times N$ matrix. Is that right? – Ben Grossmann Aug 07 '23 at 14:33
  • Regarding numerical rank deficiency, you might have some luck using some kind of preconditioner. – Ben Grossmann Aug 07 '23 at 14:39
  • @BenGrossmann: Yes, $A(\alpha)$ is a $N \times N$ matrix whose entries depend on the parameter $\alpha \in \mathbb{R}$. I already expected that preconditioning might be the way to go but I am super unfamiliar with it, so far. Most literature I read seems to be based on solving eigenvalue problems but not on null space computation (eigenvalue zero). Any suggestions are again highly welcomed. – Michael Steinke Aug 08 '23 at 12:06
  • I suspect that I'm not the first person to think of this, but I think that the following should work: if there is some $\alpha_1$ near your problematic $\alpha_0$ such that $A(\alpha_1)$ is at least numerically invertible, then $[A(\alpha_1)]^{-1}A(\alpha_0)$ should be relatively well-conditioned. The idea is that the continuity of $\alpha \mapsto [A(\alpha_1)]^{-1}A(\alpha)$ ensures that the matrix will be "close" to the identity matrix when $\alpha$ is close to $\alpha_1$. – Ben Grossmann Aug 08 '23 at 12:59
  • @BenGrossmann I'll check that and let you know. Thanks! – Michael Steinke Aug 08 '23 at 13:04

1 Answers1

1

I'm not an expert, but worked on something similar recently, so will give my 2 cents.

If rank(A) = N, then you will have a unique solution, the trivial one. If rank(A) < N, then you have an infinite number of solutions.

So you want A to be rank-deficient at $\alpha_0$, if you want a non-trivial solution.

You can calculate a matrix Z using scipy.linalg.null_space(). It has size [N, N-rank(A)].

All $\vec{x} = Z\vec{u}$ are solutions to your equation, where $\vec{u}$ is a vector of length N-rank(A).

The fact that the lowest singular values drop below the threshold is good, it's actually how the null_space() function determines what vectors from V to select. And thereby the size of the null-space. You can look at the source-code

So, to find $\alpha_0$, perhaps you can scan over different values of $\alpha$ and see when null_space(A) is not empty.

Alternatively, scan over different values of $\alpha$, do the SVD on each step and see when the lowest singular value $\sigma$ falls below some threshold(Can use the same definition as in the null_space() function). Good thing about this approach is that $\sigma(\alpha)$ is differentiable

Hope that helps. Sorry if I misunderstood xD.

Aleksk89
  • 51
  • 3