safe_inv should have a different tol from pullback_gaugetol

[Confusio]

When applying safe_inv to a vector A with a very small maximum entry, the tolerance should be scaled by the vector’s infinity norm (i.e. set to tol * norm(A, Inf)) rather than using the default pullback_gaugetol value. Otherwise, you would need to renormalize A prior to calling safe_inv.

[lkdvos]

Could you elaborate a bit on what you mean? the pullback_gaugetol implementation does exactly what you describe, that tolerance is already scaled by the norm?

[Confusio]

When the norm n is small, the pullback_gaugetol is set as eps(eltype(n))^(3 / 4) since max(n, one(n))>=1 in eps(eltype(n))^(3 / 4) * max(n, one(n)). While safe_inv need eps(eltype(n))^(3 / 4) * n, pullback_gaugetol need eps(eltype(n))^(3 / 4) * max(n, one(n)).

[lkdvos]

Ah I see what you mean now, thanks for clarifying.
I am also okay with simply setting the pullback_gaugetol to eps(n)^(3/4), that would also solve the issue I think.
@Jutho do you know if there was a particular reason to clamp the tol like this?

[Confusio]

Maybe I know the reason for this. As in TensorKit, I report the warning issue for svd AD with a small norm where the tol was set as eps(eltype(n))^(3 / 4) * n. If pullback_gaugetol is scaled by the norm, the warning will appear every time. That was the reason for pullback_gaugetol.

[Confusio]

Ah I see what you mean now, thanks for clarifying. I am also okay with simply setting the pullback_gaugetol to eps(n)^(3/4), that would also solve the issue I think.

When the norm is large, eps(n)^(3/4) will be too small. I think the current setting for pullback_gaugetol is reasonable. Maybe we need a different setting of tol for safe_inv.

[Jutho]

The pullback rules here are indeed just copied from TensorKit (and are currently not yet being used, nor is anything in this package at this point).

As for the default tolerance, its value has been trial and error so far.

[Jutho]

I am now working at the chainrules in #11 , and have been adding different keyword arguments for the different tolerances. @Confusio , maybe you can take a look if that would be more in line with what you had in mind?

[Confusio]

I suggest the following default

tol::Real=eps(eltype(S))^(3 / 4)*norm(S),
rank_atol::Real = tol,
degeneracy_atol::Real=tol
gauge_atol::Real=eps(eltype(S))^(3 / 4)*max(norm(S),one(norm(S)))

The setting tol::Real=default_pullback_gaugetol(S), where (default_pullback_gaugetol(S)=eps(eltype(S))^(3 / 4)*max(norm(S),one(norm(S)))), has problem when S has a norm far away from 1, as the safe_inv is not safe anymore and mask can't contain the expected mask indices.

[lkdvos]

Is there a reason to include the one in there? I guess that will always fail the pathological examples where you simply have a matrix which is A = 1e-15 * rand(n, n), which I think should otherwise still work. Basically, can we get away with eps(norm(S))^(3 / 4)?

Alternatively, we can also do something like isapprox:

isapprox(x, y; atol::Real=0, rtol::Real=atol>0 ? 0 : √eps, ...)

ie by default we only use a relative tolerance, but allow for some way of providing both atol and rtol?

[Confusio]

one is for the gauge gauge_atol as gauging part is aUΔU = rmul!(UΔU - UΔU', 1 / 2) which should be compared with 1. For the masking and save_inv, we need a relative tolerance. For the gauge warning gauge_atol, we need one for small norm and rescale by the norm for a large norm.

[Jutho]

I've currently merged the AD rule PR with the old defaults, as I forgot to look at this. However, I will study the argumentation made here and then incorporate this as well.

[Confusio]

Another suggestion is to use the broadening method to estimate the Δgauge and compare it with ΔA.

Δgauge = norm(
    (view(aUΔU, mask) + view(aVΔV, mask)) .* view(Sr' .- Sr, mask) ./
    (view(Sr' .- Sr, mask) .^ 2 .+ degeneracy_atol),
)

Δgauge < norm(ΔA) * gauge_atol ||
    @warn "`svd` cotangents sensitive to gauge choice: (|Δgauge| = $Δgauge)"

In general gauge_atol=1.0-4 is enough in practice for optimization.

[Jutho]

Yes, that might be a good suggestion. It is clear that rank_atol and degeneracy_atol refer to the properties of S and should thus be relative to the norm of S. gauge_atol on the other hand refers to properties of the cotangents, and should either be relative to the size thereof, or indeed use something like your suggestion. The only "problem" with the Δgauge < norm(ΔA) * gauge_atol suggestion is that this is not an absolute but relative tolerance, so maybe I should just switch everything to relative tolerances? In this case, we cannot compute the corresponding absolute tolerance, since ΔA is only known at the end of the computation.

[Confusio]

We can move the warning at the end.

[Jutho]

Yes I realize that the warning can be printed after ΔA is computed; it is more a semantic issue that your use of gauge_atol is not an absolute tolerance (which is what the a in atol is for, but a relative tolerance). So maybe switching to relative tolerance everywhere is the way to go, or use both indeed.