SVD

Stream: maths

[Hans Parshall]

@Daniel Packer and I see a route to formalizing the singular value decomposition (SVD) by following Axler's book a bit further. We would build on the work by @Frédéric Dupuis and @Heather Macbeth on adjoints and spectral theory, and then we'd like to extract SVD for real/complex matrices. Is anyone else working on this?

I should add that Dan and I are Lean novices, and so our "first draft" will probably need some help to make it useful for mathlib.

[Frédéric Dupuis]

Sounds like a great project! I'll take the opportunity to bring up again a relevant issue that came up recently: it's currently very unpleasant to work with the standard basis (i.e. pi.basis_fun, and related things like matrix.to_lin') on pi_Lp and euclidean_space. The reason is that pi.basis_fun is defined as the standard basis on ι → 𝕜 and while those other two spaces are type copies of this, Lean never figures it out, so we always have to feed in the types manually.

One can think of a few solutions, for example we could duplicate the API of pi.basis_fun, matrix.to_lin' and friends to the type copies, but I'm afraid of the code duplication this would entail. Another possibility would be to have a [has_standard_basis ι 𝕜 E] typeclass, which sounds more promising to me.

Any thoughts?

[Yaël Dillies]

Did you think about having the two "identity" equivalences as is done with docs#order_dual.to_dual / docs#order_dual.of_dual and docs#to_lex / docs#of_lex?

[Yaël Dillies]

Works great in order theory, once you fight pattern-matching a bit (cf file#data/sum/interval.lean) when dealing with inductives.

[Heather Macbeth]

I actually have been thinking of a more radical proposal, which is to demote the current metric space instance on pi (with the sup-norm) to a locale. This would allow for a rival locale, used in the inner product space library, in which it's the l2-norm.

[Frédéric Dupuis]

Yaël Dillies said:

Did you think about having the two "identity" equivalences as is done with docs#order_dual.to_dual / docs#order_dual.of_dual and docs#to_lex / docs#of_lex?

That could be another option, though I suspect having to use those equivalences all the time could get old pretty fast.

[Frédéric Dupuis]

Heather Macbeth said:

I actually have been thinking of a more radical proposal, which is to demote the current metric space instance on pi (with the sup-norm) to a locale. This would allow for a rival locale, used in the inner product space library, in which it's the l2-norm.

I'm not sure I'm getting this right -- you would propose to get rid of euclidean_space and pi_Lp entirely and replace them by locales?

[Heather Macbeth]

I would keep the type synonym pi_Lp since (at this stage) I don't see the other Lp spaces being used often enough to give headaches. But yes, I'd get rid of euclidean_space. We'd instead have that in the locale inner_product_space the metric and norm instances on finitely-indexed pi-types, including ι → 𝕜, be given by the L2 norm.

[Heather Macbeth]

It is indeed a bit of a drastic proposal, and I haven't tried it out yet.

[Frédéric Dupuis]

If we're going to fix this, I think we should also fix it for pi_Lp, it'll be used with bases sooner or later. If we go down the route of duplicating the API, we could get rid of euclidean_space, put an inner product space structure on pi_Lp 2, and hope we don't have to worry about a third type copy down the road.

[Heather Macbeth]

Frédéric Dupuis said:

with bases

What do you mean by this, sorry?

[Frédéric Dupuis]

I mean we'll want to use matrix.to_lin' and friends with those spaces at some point.

[Heather Macbeth]

Do you think so? I never have ...

[Frédéric Dupuis]

Sure, I use the 1-norm all the time and matrices are involved.

[Heather Macbeth]

Frédéric Dupuis said:

If we go down the route of duplicating the API, we could get rid of euclidean_space, put an inner product space structure on pi_Lp 2, and hope we don't have to worry about a third type copy down the road.

And with this proposal, I don't quite follow? We already have an inner product space structure on pi_Lp 2, docs#pi_Lp.inner_product_space

[Frédéric Dupuis]

Hmm, so why do we have euclidean_space again?

[Heather Macbeth]

pi_Lp is for a finite dependent family E : ι → Type*, euclidean_space is its special case E = λ i, 𝕜.

[Heather Macbeth]

Frédéric Dupuis said:

Sure, I use the 1-norm all the time and matrices are involved.

The 1-norm is also somewhat special, though (maybe we would have a third locale for it!), the other $p$-norms I personally only tend to use in some kind of auxiliary step building up to a better norm.

[Frédéric Dupuis]

Then making euclidean_space notation instead of a type copy would be a good first step!

[Heather Macbeth]

I don't think it would fix problems like in #11551, though ... a type synonym of a type synonym is still as annoying as a type synonym, I think.

[Heather Macbeth]

And there are also a lot of annoying unification errors that tend to strike when you work with a dependent pi-type, it may be that euclidean_space was introduced to circumvent that issue.

[Heather Macbeth]

Still, it would be interesting to try it!

[Heather Macbeth]

I guess with my other proposal, my concern is that I don't see how to have a "$p$-dependent locale"

[Frédéric Dupuis]

Heather Macbeth said:

I don't think it would fix problems like in #11551, though ... a type synonym of a type synonym is still as annoying as a type synonym, I think.

It still means that we could get away with only two versions of matrix.to_lin' and friends. At 3 copies I think it starts being excessive.

[Frédéric Dupuis]

In any case, I still think the [has_standard_basis ι 𝕜 E] class looks like the best way out. Is there an obvious drawback?

[Heather Macbeth]

Can you explain this idea a bit more? I didn't understand when you mentioned it before.

[Frédéric Dupuis]

We would define, say, matrix.to_lin' so that it constructs a linear map on E if the class [has_standard_basis ι 𝕜 E] is there.

[Heather Macbeth]

It's a data class, containing precisely an ι-indexed basis for E?

[Frédéric Dupuis]

Yeah

[Heather Macbeth]

I think this also seems like an elegant solution, yes. I'd be curious to go a little way with both this and the separate-locales method, and see how they work.

[Heather Macbeth]

I wonder whether with has_standard_basis you'd lose convenient things like the ability to directly consider a vector constructed using the ![3, 2, 5] notation as an element of the space.

[Heather Macbeth]

Or, if you redefined the ![3, 2, 5] notation to apply to has_standard_basis, whether you'd lose its computational power for the original setting of ι → 𝕜.

[Heather Macbeth]

We'd have to try it to find out!

[Frédéric Dupuis]

You would if you directly work with euclidean_space or ι → 𝕜. But admittedly we would end up doing this less often if we go down that route.

[Heather Macbeth]

Hmm, so what would we do with a file like
https://leanprover-community.github.io/mathlib_docs/number_theory/modular.html
which is full of explicit calculations about 2x2 matrices?

[Frédéric Dupuis]

Are we converting back and forth between matrices and linear maps there? At first glance it looks mostly like we're considering the space of matrices as a vector space.

[Heather Macbeth]

Mostly it's individual matrices, so simpler really than either of those things. But the ![3, 2, 5] notation is used quite a bit.

[Frédéric Dupuis]

But this wouldn't affect the ability to define matrices.

[Frédéric Dupuis]

Actually, could [has_coe (ι → 𝕜) E] work instead of creating a new typeclass?

[Heather Macbeth]

I see, so you would go for this
Heather Macbeth said:

I wonder whether with has_standard_basis you'd lose convenient things like the ability to directly consider a vector constructed using the ![3, 2, 5] notation as an element of the space.

(with the vector/matrix notation unaffected, but not applying to euclidean_space) rather than this
Heather Macbeth said:

Or, if you redefined the ![3, 2, 5] notation to apply to has_standard_basis, whether you'd lose its computational power for the original setting of ι → 𝕜.

[Frédéric Dupuis]

I think we should test this has_coe thing first, if it works then I think it solves all our problems.

[Heather Macbeth]

Frédéric Dupuis said:

Actually, could [has_coe (ι → 𝕜) E] work instead of creating a new typeclass?

Hmm ... I think it works to allow you to easily consider ![3, 2, 5] as an element of E, but not to allow you to easily consider ![[3, 2, 5], [0, 0, 1], [1, 1, 1]] as an operator on E?

[Frédéric Dupuis]

Well, this notation would just work like does now.

[Heather Macbeth]

It's a bit tricky, too, right? Because you don't want just any coercion, you want a coercion which is an isomorphism of vector spaces.

[Heather Macbeth]

Frédéric Dupuis said:

Well, this notation would just work like does now.

But the situation now is unsatisfactory, that's why we have all these problems with matrix.to_lin' on euclidean_space.

[Frédéric Dupuis]

Oh, does that notation not work with euclidean_space?

[Heather Macbeth]

I believe not. The notation makes a matrix (i.e., ι → ι → 𝕜), matrix.to_lin' turns a matrix into an operator on ι → 𝕜, and as we saw in #11551 it's quite painful to identify that with an operator on euclidean_space.

[Frédéric Dupuis]

Heather Macbeth said:

It's a bit tricky, too, right? Because you don't want just any coercion, you want a coercion which is an isomorphism of vector spaces.

It might be just enough to work for the few types we care about here. It'll produce garbage for a general coercion, but we could probably live with it.

[Heather Macbeth]

It's worth a shot! I think we should do some experiments.

[Yaël Dillies]

Frédéric Dupuis said:

Yaël Dillies said:

Did you think about having the two "identity" equivalences as is done with docs#order_dual.to_dual / docs#order_dual.of_dual and docs#to_lex / docs#of_lex?

That could be another option, though I suspect having to use those equivalences all the time could get old pretty fast.

It's surprisingly not. Just a matter of inserting them at the correct places, which doesn't happen too often. And I think you should have them in all cases, because they're roughly the only way to make stuff typecheck without breaking the API barrier and ending up with nonsensical goals.

[Frédéric Dupuis]

Yaël Dillies said:

It's surprisingly not. Just a matter of inserting them at the correct places, which doesn't happen too often. And I think you should have them in all cases, because they're roughly the only way to make stuff typecheck without breaking the API barrier and ending up with nonsensical goals.

Part of the reason why I'm skeptical is that in this particular case, it involves dealing with things like:

def linear_map.to_matrix' : ((n → R) →ₗ[R] (m → R)) ≃ₗ[R] matrix m n R

on type copies of n → R and m → R. So even just figuring out how to insert those equivs would be quite painful.

[Eric Wieser]

In #18574 I add matrix.to_euclidean_lin : matrix m n 𝕜 ≃ₗ[𝕜] (euclidean_space 𝕜 n →ₗ[𝕜] euclidean_space 𝕜 m) to make this a bit less painful

[Eric Wieser]

Re-reading the thread I see that I suggested this a year ago, but I think we have more lemmas about docs#pi_Lp.equiv now than we used to, so it's less painful.