find_start_pair failing on some linearly-parametrized square systems
timduff35 opened this issue · 2 comments
The code below (based off a recently-started research project with undergrads) produces a square parametric system in unknowns X, Y
and parameters A,B,C
. It is one way of formulating the Eckart-Young problem of approximating an n x n matrix by a rank-r matrix.
When we tried solving this system with monodromy, we noticed that find_start_pair
usually produced a solution with nearly-singular Jacobian. However, we can find plenty of non-singular starting pairs using the function better_start_pair
below. This function assumes only that the system is affine-linear in the parameters. It is more-or-less the same as this Macaulay2 function.
I believe the issue here is that the particular solution of the linear system computed by find_start_pair_linear_in_params
tends to have many coordinates equal to zero. At least for this particular class of systems, these sparse parameter vectors seem to produce singular systems. The idea behind better_start_pair
is to add on a dense null-vector to this particular solution. The resulting start part seems to work as expected when passed to monodromy_solve
.
This strategy ought to work for any cases where find_start_pair
already works, and I think only a few new lines of code would need to be added. Do you see any potential downsides?
using HomotopyContinuation
using LinearAlgebra
frob(X) = sum([X[i,j]^2 for i = 1:size(X)[1], j = 1:size(X)[2]]);
function ey_system(n,r)
@assert r <= n;
@var A[1:n,1:n], B[1:r, 1:n], C[1:r, 1:r];
@var X[1:n, 1:r], Y[1:n, 1:r]; # variables
critical_equations = differentiate(frob(A - X*Y'), [vec(X); vec(Y)]);
indices_removed = [[i,j] for i=1:r, j=1:r];
indices_kept = setdiff(1:2*n*r, indices_removed);
parametric_system = System(vcat(critical_equations[indices_kept], vec(B*X - C)); parameters = vcat(vec(A), vec(B), vec(C)));
num_equations = length(expressions(parametric_system));
num_variables = length(expressions(parametric_system));
@assert(num_equations == num_variables);
parametric_system
end
function better_start_pair(F)
n = length(variables(F));
x0 = vec(randn(n,1) + im * randn(n,1));
m = length(parameters(F));
idmat = convert(Matrix{ComplexF64}, I(m));
b = F(x0,zeros(m));
A = hcat([F(x0,idmat[:,i])-b for i=1:m]...);
N = nullspace(A);
weights = randn(size(N,2),1) + im * randn(size(N,2),1);
pHom = vec(N * weights);
pPart = - A \ b;
p0 = pPart + pHom;
return (x0, p0)
end
# Example
F = ey_system(8,3);
# built-in start pair
(x0, p0) = find_start_pair(F);
minimum(svdvals(jacobian(F, x0, p0))) # gives me ~ 1e-16
# "custom" start pair
(x1, p1) = better_start_pair(F);
minimum(svdvals(jacobian(F, x1, p1))) # gives me >> 1e-16
monodromy_solve(F, x1, p1) # gives me 56 = 8C3 solutions
Which HC version are you using? We recently merged this PR #532 . Does this address the issue?
Oh, oops! This was exactly the same issue, it should work now.