Ill-conditioned solve with Inverse Geometric Hodge 1

[GeorgeR227]

Resolving the issues exposed in #177 has exposed a potential problem with the dual Laplacian. I have formulated the following MWE that demonstrates the error is likely arising in Geometric Hodge 1. This is based off a testcase in test/Operators.jl. Essentially, we want to the take the Laplacian of the scalar function 2x. This should return 0, namely on the interior since the d0 is known to have boundary issues.

As shown below the formulation using Geometric Hodge is completely off and does not seem to improve with better resolution. Using the signed Hodge 2 improves the mean result but I believe this is a fluke for two reasons. 1) The standard deviation in the answer is still quite high, indicating significant non-zero entries and 2) using the Diagonal Hodge produces the correct solution of 0. In short, the signed Hodge 2 seems to have merely masked a deeper problem apparent in the Dual Laplacian on 0-forms.

using CombinatorialSpaces
using Statistics

s = triangulated_grid(100,100,10,10,Point2d);
sd = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2d}(s);
subdivide_duals!(sd, Barycenter());
interior_tris = setdiff(triangles(sd), boundary_inds(Val{2}, sd))
dd0 = dec_dual_derivative(0, sd);
d1 = dec_differential(1,sd);
hs2 = dec_hodge_star(2, sd, GeometricHodge());

twoX = map(p -> 2*p[1], sd[sd[:tri_center], :dual_point])

# Geometric Hodge
# Mean ≈ 0.0474, Std ≈ 0.298
ihs1 = dec_inv_hodge_star(1, sd, GeometricHodge());
nil = hs2 * d1 * ihs1(dd0 * twoX)
abs(mean(nil[interior_tris]))
std(nil[interior_tris])

# With incorrect signed Hodge 2
# Mean ≈ 1.4e-16, Std ≈ 0.30
nil = sign(2, sd) .* nil
abs(mean(nil[interior_tris]))
std(nil[interior_tris])

# Diagonal Hodge
# Mean ≈ 2e-16, Std ≈ 1e-15
ihs1 = dec_inv_hodge_star(1, sd, DiagonalHodge());
nil = hs2 * d1 * ihs1 * dd0 * twoX
abs(mean(nil[interior_tris]))
std(nil[interior_tris])

The cause of the error is unknown as of now but it is possible a similar sign/orientation issue may be present.

[GeorgeR227]

I've been testing this more and it seems like the Geometric Hodge 1 itself is okay. I compared the test cases present in Ayoub et al (2020) in section 3.2 and confirmed that we matched with all of those intended results. This was also tested with a variety of edge and triangle orientations and none of these affected the final result, as intended.

Considering then that we used the inverse Geometric Hodge 1 above, I noticed that the determinant of the Hodge matrix was numerically 0.

julia> det(dec_hodge_star(1, sd, GeometricHodge()))
4.992931430538642e-90

I'm not sure if this matrix is meant to be singular but I would expect it not to be since this is run on a triangulated_grid. This result leads me to believe that this is a poorly conditioned solve which would introduce error and thus cause the issues above.

When I do a roundtrip for a simple solve, I am able to recover the full solution so this I'm not sure if this is a actual problem with our implementation or if this particular example is not well-behaved.

using CombinatorialSpaces

s = triangulated_grid(100,100,10,10,Point2d);
sd = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2d}(s);
subdivide_duals!(sd, Barycenter());

invhdg1 = dec_inv_hodge_star(1, sd, GeometricHodge())
hdg_1 = dec_hodge_star(1, sd, GeometricHodge())

# This result is as intended
all(invhdg1(hdg_1 * ones(ne(sd))) .≈ -1)