search

JuliaApproximation / DomainSets.jl

A Julia package for describing domains as continuous sets of elements
MIT License
72 stars 12 forks source link

Implicit pinv of vectors #82

Closed antoine-levitt closed 3 years ago

antoine-levitt commented 3 years ago

Hi,

Investigating https://github.com/JuliaLang/julia/pull/40758 I found out that this package fails under that PR, which undefines number / vector (which makes an implicit pseudo-inverse of the vector). I haven't been able to pinpoint where exactly it fails in this package, but it's something to do with test_generic_map(m). Can you confirm if this use of the pseudo-inverse is intended, or a bug and number ./ vector is meant?

daanhb commented 3 years ago

That is definitely not intended :-) I'll try to reproduce it and I'm looking into it, but from the sound of it the aim in DomainSets would be to use pinv explicitly

daanhb commented 3 years ago

On closer look, it's somewhat ambiguous. It happened in a code path in which I know that the linear system I'm solving is overdetermined, so I'm looking for a least squares solution.

The offending line is here. The method is computing the left inverse of a linear map (i.e. y=Ax with rectangular A), evaluated at a point x. For the left inverse of the map itself, I'm using pinv here.

However, in the offending line I'm evaluating the left inverse at a point, so I don't need the full inverse. In the same way that A\b is preferred over inv(A)*b, here I was using A\b for the least squares solve rather than pinv(A)*b. I guess I was hoping that it might be faster.

I've now removed that assumption, and in any case A\b is not faster than pinv(A)*b if \ ends up doing a pinv behind the scenes. But I'm not sure that it was an invalid thing to do. I'll give it some more thought.

antoine-levitt commented 3 years ago

Thanks for investigating! It is indeed equivalent to the pinv. It's completely valid usage, sorry for the noise!

antoine-levitt commented 3 years ago

(for context: I'm investigating whether it makes sense to remove all implicit pseudo-inverses of vectors, because it is very easy to write 1 / x instead of 1 ./ x. What I'm seeing is that many people including you use x\y for least squares, so in this context it makes sense to have implicit pseudo-inverses)

daanhb commented 3 years ago

I do sympathise with the cause :-)