Singular Value Decomposition

Stream: Is there code for X?

[Niels Voss]

Is there code for or a plan for singular value decomposition (of matrices or linear maps)? I found #6042 by @MohanadAhmed, but it seems that it's not general enough to be merged as-is. I was wondering if anything had replaced it yet, or if there was a plan to finish it. If not, what's the ideal generalization of SVD for mathlib (I assume there are ways to generalize SVD for finite-dimensional inner product spaces to maybe something that is not finite-dimensional)?

[Jireh Loreaux]

There the Schmidt decomposition for compact operators on a Hilbert space, but that's different enough that we should have both, not try to shoehorn SVD to fit that mold.

[Niels Voss]

Is the schmidt decomposition in mathlib, or are you just saying that it would be a candidate for a generalization of SVD?

[Niels Voss]

Can the schmidt decomposition be generalized further to infinite dimensions (I only saw the finite-dimensional case on wikipedia) or to larger tensor products than just two spaces, or is it as generalized as it will probably get?

[Niels Voss]

Would singular value decomposition of linear maps between finite-dimensional inner product spaces be on topic for mathlib? Or would it need to be generalized to the infinite-dimensional case? We only have spectral theorem for linear maps between finite-dimensional spaces

[Niels Voss]

I was considering something like

variable {𝕜 : Type*} [RCLike 𝕜]
variable {E : Type*} [NormedAddCommGroup E] [InnerProductSpace 𝕜 E]
variable {F : Type*} [NormedAddCommGroup F] [InnerProductSpace 𝕜 F]
variable {n : ℕ} (hn : Module.finrank 𝕜 E = n)
variable (T : E →ₗ[𝕜] F)
variable [FiniteDimensional 𝕜 E] [FiniteDimensional 𝕜 F]

noncomputable def LinearMap.singularValues (i : Fin n) : ℝ :=
  Real.sqrt ((LinearMap.isSymmetric_adjoint_comp_self T).eigenvalues hn i)

(where LinearMap.isSymmetric_adjoint_comp_self comes from #Is there code for X? > T*T is positive when domain and codomain are different) This definition has the nice property of preserving multiplicities of singular values, as opposed to something like a HasSingularValue predicate, which wouldn't.

[Yaël Dillies]

Niels Voss said:

Would singular value decomposition of linear maps between finite-dimensional inner product spaces be on topic for mathlib? Or would it need to be generalized to the infinite-dimensional case?

Would generalising to the infinite-dimensional case involve throwing away (most of) the finite-dimensional proof? If not, then mathlib will be happy to have finite-dimensional SVD.

[Niels Voss]

To be honest, I don't personally know how to generalize SVD to the infinite-dimensional case. Do you know what the "best" generalization is?

[Niels Voss]

Also, for T : E -> F, when mathematicians say "list of singular values", does this list generally

• have length rank(T), so that all the singular values are positive
• have length min(dim(E), dim(F))
• have length dim(E) (as it does in LADR 4th edition)
• have countably infinite length, but all but finitely many of the singular values are zero?

I'm trying to decide whether indexing the list of singular values by Fin n or Nat is a better idea (and if the latter, whether or not I should Finsupp)

[Niels Voss]

There's also this alternative formulation from https://en.wikipedia.org/wiki/Singular_value
image.png

[Niels Voss]

I changed it to this definition:

variable {𝕜 : Type*} [RCLike 𝕜]
variable {E : Type*} [NormedAddCommGroup E] [InnerProductSpace 𝕜 E]
variable {F : Type*} [NormedAddCommGroup F] [InnerProductSpace 𝕜 F]
variable (T : E →ₗ[𝕜] F)
variable [FiniteDimensional 𝕜 E] [FiniteDimensional 𝕜 F]

noncomputable def LinearMap.singularValues : ℕ →₀ ℝ :=
  Finsupp.embDomain Fin.valEmbedding <|
    Finsupp.ofSupportFinite
      (Real.sqrt ∘ (LinearMap.isSymmetric_adjoint_comp_self T).eigenvalues rfl)
      (Set.toFinite _)

with the characterization

theorem LinearMap.apply_singularValues_of_lt {n : ℕ} (h : n < Module.finrank 𝕜 E)
  : T.singularValues n
    = √((LinearMap.isSymmetric_adjoint_comp_self T).eigenvalues rfl ⟨n, h⟩) := sorry

theorem LinearMap.apply_singularValues_of_le {n : ℕ} (h : Module.finrank 𝕜 E ≤ n)
  : T.singularValues n = 0 := sorry

theorem LinearMap.apply_singularValues_fin {n : ℕ} (hn : Module.finrank 𝕜 E = n) (i : Fin n)
  : T.singularValues ↑i = √((LinearMap.isSymmetric_adjoint_comp_self T).eigenvalues hn i) := sorry

[Sébastien Gouëzel]

I'd rather allow infinite dimension (which means that, in general, singular values don't form a finitely supported sequence). For instance, the singular values of the Sobolev injections are interesting objects, as they quantify "how much compact" the injection is. More generally, approximation numbers (the inf of the norm of T - L for finite rank L) are very much relevant in this setting, as they make sense for Sobolev spaces not based on L^2.

[Jireh Loreaux]

I think defining approximation numbers would be great. The thing to prove, as far as I'm concerned, with respect to SVD would be that for a compact operator there are two orthonormal sequences such that the compact operator is the norm convergent sum of the rank-1 elementary tensors induced by pairing these vectors and multiplying by the approximation numbers.

[Niels Voss]

Since I don't know the actual math yet: Are approximation numbers defined for continuous linear maps only or all linear maps? And can they be infinite? Should the codomain be Real, NNReal, or ENNReal?

[Sébastien Gouëzel]

It's only defined for continuous linear maps (for non-continuous linear maps, all the approximation numbers would be infinite, so they wouldn't be interesting at all). Note that in finite dimension, all linear maps are continuous, so there's no loss in generality.

[Niels Voss]

To confirm, for non-continuous linear maps, every approximation number is infinite because subtracting a finite-rank linear map will never result it it becoming continuous, and all non-continuous maps have infinite norm, right?

[Niels Voss]

Defining it only for continuous linear maps would make it slightly more annoying to use, because the conversion between linear and continuous linear maps is not automatic even for finite dimensions. This is why we have LinearMap.adjoint in addition to ContinuousLinearMap.adjoint

[Sébastien Gouëzel]

I would argue that it's more natural to define it only for continuous linear maps, even in finite dimension: it's an analytical notion, depending on the norm, so naturally the realm of continuous linear maps.

[Jireh Loreaux]

Sébastien, I think that's only half true. People will certainly want SVD for matrices without caring much about the norm / topology.

[Niels Voss]

If the only cost is that people have to do T.toContinuousLinearMap.singularValues instead of T.singularValues it isn't that bad

[Sébastien Gouëzel]

I think that for matrices we should have a shortcut indeed. Note that matrices are not the same as linear maps: to me, one should have them for continuous linear maps because that's the right conceptual framework, and then for matrices because it's used so often.

[Niels Voss]

And nothing for regular linear maps, right?

[Jireh Loreaux]

Great, sounds like we're on the same page then.

[Kim Morrison]

I love the plan, but I hope that if there is not an implementation plan (i.e. a person saying they are doing it in a time frame) that this would not get in the way of someone contributing an implementation for matrices, which there is demand for and has been for a while. I think it should suffice in an implementation for matrices to include a module doc link to this thread, giving a pointer to the proposed upgrade path.

[Niels Voss]

Should the codomain of ContinuousLinearMap.singularValue be Real or NNReal?

[Jireh Loreaux]

I would be tempted to define it as ℕ → ℝ≥0 as subtraction plays a lesser role in the context of singular values. Nonnegativity is more important I think.

[Niels Voss]

Does this look like a reasonable definition then?

/--
When the vector spaces are not necessarily finite, these are called the approximation numbers.
Note: we use ≤ n instead of < n because we start the singular values at 0, not 1, so if we did
< n the first singular value would be repeated.
-/
noncomputable def ContinuousLinearMap.singularValue {𝕜 : Type*} [NontriviallyNormedField 𝕜]
  {X : Type*} [SeminormedAddCommGroup X] [NormedSpace 𝕜 X]
  {Y : Type*} [SeminormedAddCommGroup Y] [NormedSpace 𝕜 Y]
  (T : X →L[𝕜] Y) (n : ℕ) : ℝ≥0 :=
  -- NOTE: We can't use `⨅ H ∈ {...}` because of
  -- https://leanprover-community.github.io/extras/pitfalls.html#accidental-double-iinf-or-isup
  ⨅ H : {S : X →L[𝕜] Y // S.rank ≤ n}, ‖T - H‖₊

[Niels Voss]

What way would you recommend to define singular values for Matricies? For a matrix M, I can either say that

• The singular values of M are the singular values of the operator T produced by converting M to a linear map between EuclideanSpace k m and EuclideanSpace k n (probably with docs#Matrix.toEuclideanLin). This basically the approach Matrix.IsHermitian.eigenvalues₀ uses. Requires M to be finite (or at least needs M to have a finite number of columns) and the field to be RCLike.
• The singular values of M are defined by inf {||M - N|| : rank(N) < n} where we use the norm defined in Mathlib.Analysis.CStarAlgebra.Matrix. Requires M to be finite and the field to be RCLike.
• The singular values of M are defined by the square roots of the eigenvalues of T*T, produced by docs#Matrix.IsHermitian.eigenvalues₀. Requires M to be finite and the field to be RCLike.

All of these approaches only work for finite-dimensional Real and Complex. Between these three approaches, which is the best, and is there a different approach that can be used? If only finite matrices are supported, should the definition be in Finsupp, like above?

[Niels Voss]

The first option would probably be something like

open NNReal in
noncomputable def Matrix.singularValue {𝕜 : Type*} [RCLike 𝕜]
  {n m : Type*} [Fintype n] [Fintype m] [DecidableEq m] (M : Matrix n m 𝕜)
  : ℕ → ℝ≥0 := M.toEuclideanLin.toContinuousLinearMap.singularValue

or

open NNReal in
noncomputable def Matrix.singularValue {𝕜 : Type*} [RCLike 𝕜]
  {n m : Type*} [Fintype n] [Fintype m] [DecidableEq m] (M : Matrix n m 𝕜)
  : ℕ →₀ ℝ≥0 :=
  Finsupp.ofSupportFinite
    M.toEuclideanLin.toContinuousLinearMap.singularValue
    sorry

(I can probably get rid of the [DecidableEq m] since it is noncomputable anyways)

[Niels Voss]

Jireh Loreaux said:

I think defining approximation numbers would be great. The thing to prove, as far as I'm concerned, with respect to SVD would be that for a compact operator there are two orthonormal sequences such that the compact operator is the norm convergent sum of the rank-1 elementary tensors induced by pairing these vectors and multiplying by the approximation numbers.

Sorry if this is a stupid question, but since I don't know the field that well, can you confirm whether this holds true between arbitrary normed spaces, arbitrary Banach spaces, Banach spaces that satisfy the Approximation Property, or just Hilbert spaces? And would a proof for just hilbert spaces be accepted, or would it need to be as general as possible? My understanding is that the standard proof for Hilbert spaces depends on the spectral theorem for Hilbert spaces, which does not yet exist in mathlib. Also, does this sum converge conditionally, unconditionally but not absolutely, or absolutely?

[Jireh Loreaux]

The fact about the existence of orthonormal sequences only holds in the context of Hilbert spaces. @Frédéric Dupuis and @Heather Macbeth at one point had the spectral theorem for compact (selfadjoint? probably not normal) operators, but I think it may be Lean 3 code and porting it was never finished. It would be nice to make that available. The sum converges in operator norm.

[Frédéric Dupuis]

Yes, I haven't forgotten about this! :-) Indeed it's Lean 3 code, and it was written before we had compact operators in mathlib, so it's based on a less general definition that had just enough API to get through this proof. Every once in a while I tell myself I should port it, and then I get distracted by other things.

[Niels Voss]

In that case, is SVD most likely blocked behind the spectral theorem? I don't mind postulating the spectral theorem for the SVD so that they can be worked on in parallel.

[Niels Voss]

I have the following extremely rough draft of the statement of SVD for matrices (assuming singular values have already been defined for linear maps as above). The names and which parameters are implicit will certainly need to change. But I want to know if I'm at all on the right track or not.

The idea is that if we are clever, we can make the standard SVD, thin SVD, compact SVD, and maybe even truncated SVD all the same definition. Given a m x n matrix M, we decompose M into M = UΣV* where

• U is m x a
• Σ is a x b
• V* is b x n

for some a and b.

• If a = m and b = n, then this gives the standard SVD.
• If a = b = min(m,n), then this gives the thin SVD.
• If a = b = rank(M), then this gives the compact SVD.
• If a = b < rank(M), then this gives the truncated SVD.

However, I wasn't sure how to handle the fact that m,n,a,b are arbitrary, not-necessarily ordered, finite types. To do this, I postulate injective functions f : a -> Nat and g : b -> Nat, which will normally be set to Fin.val in case a and b are Fin _. An alternative option would be to postulate r : Nat and f : Fin r -> a, g : Fin r -> b.

In any case, I would like to know if this is at all a reasonable approach or if something major needs to be changed.

open NNReal in
noncomputable def Matrix.singularValue {𝕜 : Type*} [RCLike 𝕜]
  {n m : Type*} [Fintype n] [Fintype m] [DecidableEq m] (M : Matrix n m 𝕜)
  : ℕ → ℝ≥0 :=
  M.toEuclideanLin.toContinuousLinearMap.singularValue

-- TODO: Fix names
noncomputable def Matrix.leftsing {𝕜 : Type*} [RCLike 𝕜]
  {m n : Type*} [Fintype m] [Fintype n] [DecidableEq n]
  (a : Type*) [Fintype a] (f : a → ℕ)
  (M : Matrix m n 𝕜)
  : Matrix m a 𝕜 := sorry

noncomputable def Matrix.rightsing {𝕜 : Type*} [RCLike 𝕜]
  {m n : Type*} [Fintype m] [Fintype n] [DecidableEq n]
  (b : Type*) [Fintype b] (g : b → ℕ)
  (M : Matrix m n 𝕜)
  : Matrix n b 𝕜 := sorry

noncomputable def Matrix.centersing {𝕜 : Type*} [RCLike 𝕜]
  {m n : Type*} [Fintype m] [Fintype n] [DecidableEq n]
  (M : Matrix m n 𝕜)
  (a : Type*) [Fintype a] (f : a → ℕ) (b : Type*) (g : b → ℕ)
  : Matrix a b 𝕜 := Matrix.of
    fun i j ↦ if f i = g j then RCLike.ofReal (M.singularValue (f i)).toReal else 0

theorem Matrix.SVD {𝕜 : Type*} [RCLike 𝕜]
  {m n : Type*} [Fintype m] [Fintype n] [DecidableEq n]
  (a : Type*) [Fintype a] (f : a → ℕ) (hf : f.Injective)
  (b : Type*) [Fintype b] (g : b → ℕ) (hg : g.Injective)
  (M : Matrix m n 𝕜)
  (hf₂ : M.singularValue.support ⊆ Set.range f) (hg₂ : M.singularValue.support ⊆ Set.range g)
  : M.leftsing a f * M.centersing a f b g * (M.rightsing b g).conjTranspose = M := sorry

[Niels Voss]

Also, most of the time when SVD is used in proofs, it seems to go like "Let UΣV* be a singular value decomposition, ...". Are there any proofs which go "Notice that UΣV* meets the criteria for a singular value decomposition of M, so it must satisfy property X"?

[Niels Voss]

@Frédéric Dupuis Is this the most recent version of your spectral theorem code? https://github.com/leanprover-community/mathlib3/tree/dupuisf/linear_map_spectral I just want to get a general sense of what the API will look like for the spectral theorem

[Frédéric Dupuis]

I think the latest one is https://github.com/leanprover-community/mathlib3/tree/semilinear-paper-update

[Niels Voss]

I've had a lot of difficulty finding what I consider "good" definitions for singular values and SVD. I think I've settled on a few definitions that I consider reasonably optimal, but would it be possible to request a review of these definitions and theorem statements before I commit a bunch of time to trying to prove them? Or would it be best to try proving them first, and then get review of the definitions during the PR review?

[Ruben Van de Velde]

You free to ask, but I suspect one of the questions will be "how easy is it to prove the first few dozen basic lemmas about the definition"

[Li Jiale]

(deleted)

[Niels Voss]

There is a more general version of SVD which is defined for arbitrary bounded operators between Hilbert spaces, as described in https://www.mdpi.com/2227-7390/8/8/1346. However, I don't think that mathlib has the necessary tools to prove this (in particular, we need spectral theorem for bounded self-adjoint operators between Hilbert spaces). Is it fine if I just work on the compact operator SVD for Hilbert spaces instead of trying to do the most general version?

[Jireh Loreaux]

1. I think above I mentioned explicitly the compact case, so yes, go with that.
2. The SVE mentioned in that paper I think is actually not particularly useful, which I argue is essentially because H →L[ℂ] H is type I as a von neumann algebra, and so what we really want is atomic (i.e., a sequence), not some measurable function.
3. In the setting of type II von Neumann algebras there is a function μ : Set.Ici 0 → ℝ≥0∞ (which is finite if the von Neumann algebra is finite) which plays the same role as the singular values. In the case when the algebra is infinite, μ (just like the singular values) is only interesting in the case when the operator is "compact" (i.e., is a norm limit of operators whose ranges have finite trace).

[Jireh Loreaux]

@Niels Voss where are your definitions exactly?

[Niels Voss]

I can post the formal definitions in a second but in the meantime here's an informal sketch that I wrote up in obsidian
image.png

[Jireh Loreaux]

seems reasonable.

[Niels Voss]

Here's what I have right now for formal definitions. (I'm not sure about the isometries in the "SVD" section, the idea is that maybe they could be used to prove a more familiar version of SVD. I'm considering deleting the ContinuousLinearMap.___SingularMatrix definitions because I think they are redundant with the Matrix.___SingularMatrix definitions.)

import Mathlib

open NNReal in
noncomputable def ContinuousLinearMap.singularValue {𝕜 : Type*} [NontriviallyNormedField 𝕜]
  {X : Type*} [SeminormedAddCommGroup X] [NormedSpace 𝕜 X]
  {Y : Type*} [SeminormedAddCommGroup Y] [NormedSpace 𝕜 Y]
  (T : X →L[𝕜] Y) (n : ℕ) : ℝ≥0 :=
  -- NOTE: We can't use `⨅ H ∈ {...}` because of
  -- https://leanprover-community.github.io/extras/pitfalls.html#accidental-double-iinf-or-isup
  -- Also, this involves a coercion Nat → Cardinal, so I might redo that part of the definition
  -- later.
  ⨅ H : {S : X →L[𝕜] Y // S.rank ≤ n}, ‖T - H‖₊

open NNReal in
noncomputable def Matrix.singularValue {𝕜 : Type*} [RCLike 𝕜]
  {n m : Type*} [Fintype n] [Fintype m] [DecidableEq m] (M : Matrix n m 𝕜)
  : ℕ → ℝ≥0 :=
  M.toEuclideanLin.toContinuousLinearMap.singularValue

section SVD

-- TODO: When I wrote up this section for the first time, I forgot that this only works for compact
-- linear maps, and not just continuous linear maps. So the following definitions will need an
-- extra hypothesis stating T is compact, and will probably live in the IsCompactOperator namespace
-- instead of the ContinuousLinearMap namespace.

variable {𝕜 : Type*} [RCLike 𝕜]
  {E : Type*} [NormedAddCommGroup E] [InnerProductSpace 𝕜 E] [CompleteSpace E]
  {F : Type*} [NormedAddCommGroup F] [InnerProductSpace 𝕜 F] [CompleteSpace F]
  (T : E →L[𝕜] F)

noncomputable def ContinuousLinearMap.rightSingularVector : ℕ → E :=
  have := T -- temporary trick to make sure that the definition includes T as a parameter
  sorry

noncomputable def ContinuousLinearMap.leftSingularVector : ℕ → F :=
  have := T
  sorry

-- TODO: Rename
open InnerProductSpace in
theorem ContinuousLinearMap.SVD (v : E)
  : T v = ∑' n : ℕ, (RCLike.ofReal ↑(T.singularValue n) : 𝕜) •
      ⟪v, T.rightSingularVector n⟫_𝕜 • T.leftSingularVector n := sorry

/--
The first `dim(E)` vectors of `ContinuousLinearMap.rightSingularVector` but bundled into an
orthonormal basis.
-/
noncomputable def ContinuousLinearMap.rightSingularOrthonormalBasis
  [FiniteDimensional 𝕜 E] {n : ℕ} (hn : Module.finrank 𝕜 E = n) : OrthonormalBasis (Fin n) 𝕜 E :=
  OrthonormalBasis.mk (v := T.rightSingularVector ∘ Fin.val) sorry sorry

/--
The first `dim(F)` vectors of `ContinuousLinearMap.leftSingularVector` but bundled into an
orthonormal basis.
-/
noncomputable def ContinuousLinearMap.leftSingularOrthonormalBasis
  [FiniteDimensional 𝕜 F] {m : ℕ} (hn : Module.finrank 𝕜 F = m) : OrthonormalBasis (Fin m) 𝕜 F :=
  OrthonormalBasis.mk (v := T.leftSingularVector ∘ Fin.val) sorry sorry

end SVD

section MatrixSVD

variable {𝕜 : Type*} [RCLike 𝕜]
  {m n : Type*} [Fintype m] [Fintype n] [DecidableEq n]
  -- a and b should be thought of as subtypes of ℕ.
  {a b : Type*} [Fintype a] [Fintype b]
  (M : Matrix m n 𝕜)

variable
  {E : Type*} [NormedAddCommGroup E] [InnerProductSpace 𝕜 E] [CompleteSpace E]
  {F : Type*} [NormedAddCommGroup F] [InnerProductSpace 𝕜 F] [CompleteSpace F]
  [FiniteDimensional 𝕜 E] [FiniteDimensional 𝕜 F]
  (T : E →L[𝕜] F)

noncomputable def ContinuousLinearMap.rightSingularMatrix
  (rb : Module.Basis n 𝕜 E) (g : b ↪ ℕ)
  : Matrix n b 𝕜 := Matrix.of fun i j ↦ rb.repr (T.rightSingularVector (g j)) i

noncomputable def ContinuousLinearMap.leftSingularMatrix
  (lb : Module.Basis m 𝕜 F) (f : a ↪ ℕ)
  : Matrix m a 𝕜 := Matrix.of fun i j ↦ lb.repr (T.leftSingularVector (f j)) i

noncomputable def ContinuousLinearMap.centerSingularMatrix
  (f : a ↪ ℕ) (g : b ↪ ℕ)
  : Matrix a b 𝕜 := Matrix.of fun i j ↦
    if f i = g j then RCLike.ofReal ↑(T.singularValue (f i)) else 0

theorem ContinuousLinearMap.matrixSVD
  (lb : Module.Basis m 𝕜 F) (rb : Module.Basis n 𝕜 E) (f : a ↪ ℕ) (g : b ↪ ℕ)
  (hf : T.singularValue.support ⊆ Set.range f) (hg : T.singularValue.support ⊆ Set.range g)
  : T.toMatrix rb lb
    = T.leftSingularMatrix lb f * T.centerSingularMatrix f g *
      (T.rightSingularMatrix rb g).conjTranspose := sorry

noncomputable def Matrix.rightSingularMatrix
  (g : b ↪ ℕ)
  : Matrix n b 𝕜 :=
  M.toEuclideanLin.toContinuousLinearMap.rightSingularMatrix (Pi.basisFun 𝕜 n) g

noncomputable def Matrix.leftSingularMatrix
  (f : a ↪ ℕ)
  : Matrix m a 𝕜 :=
  M.toEuclideanLin.toContinuousLinearMap.leftSingularMatrix (Pi.basisFun 𝕜 m) f

noncomputable def Matrix.centerSingularMatrix
  (f : a ↪ ℕ) (g : b ↪ ℕ)
  : Matrix a b 𝕜 :=
  M.toEuclideanLin.toContinuousLinearMap.centerSingularMatrix f g

theorem Matrix.SVD
  (f : a ↪ ℕ) (g : b ↪ ℕ)
  (hf : M.singularValue.support ⊆ Set.range f) (hg : M.singularValue.support ⊆ Set.range g)
  : M = M.leftSingularMatrix f * M.centerSingularMatrix f g *
    (M.rightSingularMatrix g).conjTranspose := sorry

end MatrixSVD

[Niels Voss]

The motivation for the nonstandard definition of the matrix SVD is that it is actually general enough to encompass the standard svd, thin svd, compact svd, and truncated svd at the same time.

[Niels Voss]

Here's the way it works:

Let M be an m×n matrix, and let a and b be natural numbers. Then, if rank(M) <= a and rank(M) <= b, then M = UΣV*, where:

• U is m×a (columns of U come from the left singular vectors f_1,f_2,...)
• Σ is a×b (entries of Σ come from the singular values of M)
• V* is b×n (columns of V come from the right singular vectors e_1,e_2,...)

We can recover some of the original forms of SVD as follows:

• Full SVD: Set a = m, b = n
• Thin SVD: Set a = b = min(m,n)
• Compact SVD: Set a = b = rank(M)
• Truncated SVD: Set a = b = t < rank(M), except M = UΣV* won’t hold.

[Niels Voss]

Since mathlib seems to prefer smaller PRs, could we for instance PR just the definition of the singular values and basic theorems, or would the whole SVD need to be PR'd at once?

[Jireh Loreaux]

Surely the PR for singular values can be on its own.

[Niels Voss]

Should we have a theorem which states that for the finite dimensional case, the squares of the singular values are the LinearMap.IsSymmetric.eigenvalues of T*T, or should we wait until that gets generalized to infinite dimensions?

[Jireh Loreaux]

I think it's fine to add the finite dimensional version, because the infinite dimensional version will need a compactness assumption anyway.

[Arnav Mehta]

Hey my name is Arnav I am working with @Niels Voss and we are wondering if we should get in singular values for matrices at the same time as we PR for singular values for linear maps, and also where in Mathlib folder structure should this file be placed

[Jireh Loreaux]

I would request that the singular values for matrices comes as a separate PR, especially to reduce the review burden on the initial PR. Getting definitions right is always hard, so making that initial PR do as few things as possible is better.

[Jireh Loreaux]

As for where the file should be placed, I think Analysis.Normed.Operator is probably a good folder for this.

[Niels Voss]

Sorry for the delays in working on this. Just a few super minor (though hopefully not bike-shedding) questions about the definition of singular values because these are slightly more annoying to change later:

• Should the file be SingularValues.lean or SingularValue.lean or ApproximationNumbers.lean?
• Should it be called ContinuousLinearMap.singularValues or ContinuousLinearMap.singularValue or ContinuousLinearMap.approximationNumbers?
• It currently uses the SeminormedAddCommGroup typeclass. Should it use NormedAddCommGroup instead? (I chose to use the most general one, because that seems like the Mathlib convention)
• Right now it is defined as

public noncomputable def ContinuousLinearMap.singularValue {𝕜 : Type*} [NontriviallyNormedField 𝕜]
  {X : Type*} [SeminormedAddCommGroup X] [NormedSpace 𝕜 X]
  {Y : Type*} [SeminormedAddCommGroup Y] [NormedSpace 𝕜 Y]
  (T : X →L[𝕜] Y) (n : ℕ) : ℝ≥0 :=
  ⨅ H : {S : X →L[𝕜] Y // S.rank ≤ ↑n}, ‖T - H‖₊

where the ↑n casts n to a Cardinal. But in this topic: #mathlib4 > How to express a linear map has finite rank? some people criticized the use of Cardinal. Should I attempt to avoid Cardinal in the definition?

• Regarding the module system, should the definition be exposed or not? If not, is there an attribute that automatically generates equalities?

[Kevin Buzzard]

Because the rank can be infinite, if you want to avoid cardinals you'd have to say "S has finite rank and the finrank <= n" -- what you have looks more natural to me (and less of a mouthful). You might want to make some API to avoid having to think about cardinals e.g. a theorem saying "if S has finite rank and finrank <= n then singularValue <= ||T-H||_+" or maybe even the theorem saying that it's the inf of those things.

[Jireh Loreaux]

• Analysis/Normed/Operator/SingularValue.lean seems nice.
• I think singularValue is nice for discoverability. I've only heard approximation numbers in relation to operators on Banach spaces, not for matrices.
• Seminormed is the right choice.
• Cardinal is fine because we have no notion of erank at this point.

[Niels Voss]

And singularValue is preferred over singularValues? I know that there's something like LinearMap.IsSymmetric.eigenvalues, which is plural

[Jireh Loreaux]

That's a good point. I don't have a strong opinion about pluralization.

[Levi G]

Hello everyone. My name is Levi, a Statistics Student in 3rd year going to 4th year. Previously in another PR, I had attempted to formalize SVD with the help of AI but my code was not general enough and the PR did not meet the Mathlib quality standards.

[Levi G]

However, I have made significant progress. I believe that I have something significantly much better than I had before and believe I have addressed almost most of the main issues raised by @Jireh Loreaux in a different thread.

[Levi G]

According to what I was told by @Jireh Loreaux, @Niels Voss was the one working on SVD and according to some feedback I got from others, the ideal way of working in contributing to Mathlib is by discussing design of the code with others and coming up with small snippets of code.

[Levi G]

I would like to share that progress and get your feedback especially on whether it is the right way of doing things.

[Levi G]

The code is quite large which I had been told is not the ideal way of doing things. Anyways, here it is:
T17.1 - PR Prep SVD.lean

[Niels Voss]

Hi, for the reasons discussed in this thread, the best formalization of singular values are the approximation numbers defined in #32126. My current plan is to connect these approximation numbers with the square roots of eigenvalues of T*T for a compact operator T between Hilbert Spaces, but that proof is relatively complicated, which is why the PR is currently unfinished.

[Levi G]

Niels Voss said:

Hi, for the reasons discussed in this thread, the best formalization of singular values are the approximation numbers defined in #32126. My current plan is to connect these approximation numbers with the square roots of eigenvalues of T*T for a compact operator T between Hilbert Spaces, but that proof is relatively complicated, which is why the PR is currently unfinished.

Ok, understood. Thank you for feedback. You are working on infinite dimensional case which is much more difficult. But what do you think of my proof for the finite dimensional case? (Is it at least directionally useful) I can try fix the defs especially for singular values.

[Levi G]

@Niels Voss Here is a blueprint I have started working on that is in line with your approach (though has quite a bit of sorrys). What do you think? Is it the right approach and structure?

[Levi G]

CompactOperatorSVD.lean

[Levi G]

The above:point_up: is Work-in-Progress hence the many sorrys. Would love hear your feedback.

[Johan Commelin]

What happened with your claim to drop this project?
Levi G said:

Johan Commelin Understood. It's just that I was hoping for some more feedback because this latest version is significantly better than anything that I had given before and also attempted to address their feedback. I hope you can at least take a look at this one. However, I have understood your point, it is a minefield so I will just stop pressing on this matter, hold-off completely and move on to other projects. Thank you for your feedback.

[Snir Broshi]

Another SVD implementation #mathlib4 > Matrix analysis

[Levi G]

Johan Commelin said:

What happened with your claim to drop this project?
Levi G said:

Johan Commelin Understood. It's just that I was hoping for some more feedback because this latest version is significantly better than anything that I had given before and also attempted to address their feedback. I hope you can at least take a look at this one. However, I have understood your point, it is a minefield so I will just stop pressing on this matter, hold-off completely and move on to other projects. Thank you for your feedback.

Sorry @Johan Commelin. After I coming up with the finite-dimensional SVD, I couldn't help myself. I wanted feedback on it and wanted to be of assistance to @Niels Voss who was working on the more general and difficult case of SVD for compact operators on Hilbet Spaces and its API. So I couldn't help myself, please. As of right now I have made significant progress on the general case and will soon be sharing progress if that's alright.

[Levi G]

The progress being:

All s-number axioms proven (S2, S3 left/right)

Rayleigh Quotient Theorem - Infinite dimensional case for Compact Self-Adjoint Operators

Spectral Theorem for Compact Self-Adjoint Operators on Hilbert Spaces

Compactness characterization of operators

[Levi G]

It's a big file though.

[Jireh Loreaux]

@Levi G the fact that it's a big file is part of the problem. At this point, you don't have experience contributing to Mathlib, so a contribution of this magnitude is very hard to review. If you start with smaller projects and learn good design patterns from receiving review, then it gets easier because your code is higher quality. Let me highlight one example as evidence:

theorem IsCompactOperator.eigenvalue_tendsto_zero
    (hT : IsCompactOperator T) (_hSA : IsSelfAdjoint T)
    {μ : ℕ → ℝ} {e : ℕ → E}
    (he_orth : Orthonormal 𝕜 e)
    (he_eigen : ∀ n, T (e n) = (μ n : 𝕜) • (e n))
    (hμ_antitone : Antitone (fun n ↦ |μ n|)) :
    Tendsto (fun n ↦ |μ n|) atTop (𝓝 0) := by

This theorem has nothing to do with selfadjiont operators, or real eigenvalues, or Antitone in any fundamental way. That is, what's really going on here is that since e is orthonormal, is converges weakly to zero, and compact operators take weakly convergent sequences to norm convergent ones.

Again, I would love for you to become a productive member of the community, but please start with smaller projects. It may even be the case that there are peices of your code worth salvaging, but I don't have the energy or the time to try to sift through 2000 lines of code looking for a diamond in the rough.

[Levi G]

@Jireh Loreaux Thank you for your feedback and time spent reviewing my code. Long, bad-quality code is hard to review - Understood. My new file is around 3k lines - too long.

[Levi G]

In light of the difficulty of reviewing long low quality code - I'll share what I think are the parts of interest in my new code: the Rayleigh Quotient Theorem and Spectral Theorem for Compact Self-Adjoint Operators on Hilbert Spaces. Then hold-off from the project (finally) and move on to smaller project where I can improve.

[Levi G]

/-- Key Lemma 2: For a nonzero compact self-adjoint operator, ‖T‖ or -‖T‖ is an eigenvalue.

The proof uses the Rayleigh quotient: for self-adjoint `T`, the operator norm is the
supremum of the absolute Rayleigh quotient on the unit sphere. For compact `T`, we
build a maximizing sequence on the unit sphere and use compactness to extract a
convergent subsequence, whose limit is an eigenvector.

This theorem is the compact self-adjoint case that Mathlib's `Rayleigh.lean` leaves
as a TODO; here we supply a standalone implementation using `hOpNorm_eq_sSup_absRayleigh`
(proven via the Polarization Identity) as the key spectral-theoretic input.
-/
theorem IsSelfAdjoint.exists_norm_eq_abs_eigenvalue (hSA : IsSelfAdjoint T)
    (hT : IsCompactOperator T) (hT_ne : T ≠ 0) :
    ∃ (μ : ℝ) (e : E), ‖e‖ = 1 ∧ T e = (μ : 𝕜) • e ∧ |μ| = ‖T‖ := by
  -- ═══════════════════════════════════════════════════════════════════════════
  -- PROOF STRATEGY (completing the TODO in Mathlib's Rayleigh.lean):
  -- ═══════════════════════════════════════════════════════════════════════════
  --
  -- For a nonzero compact self-adjoint operator T:
  -- 1. Either sup or inf of the Rayleigh quotient has absolute value = ‖T‖
  -- 2. Take a maximizing (or minimizing) sequence on the unit sphere
  -- 3. Use compactness to extract a convergent subsequence of T(xₙ)
  -- 4. Show the sequence xₙ itself converges using ‖Txₙ - μxₙ‖ → 0
  -- 5. The limit is an eigenvector with eigenvalue μ where |μ| = ‖T‖
  --
  -- Key insight: For self-adjoint T and maximizing sequence xₙ with ⟪Txₙ,xₙ⟫ → M:
  --   ‖Txₙ - M·xₙ‖² = ‖Txₙ‖² - 2M⟪Txₙ,xₙ⟫ + M²
  --                 = ⟪T²xₙ,xₙ⟫ - 2M⟪Txₙ,xₙ⟫ + M² → 0
  -- because ⟪T²xₙ,xₙ⟫ ≤ ‖T‖·⟪Txₙ,xₙ⟫ → ‖T‖·M and M = ‖T‖ or M = -‖T‖.
  -- ═══════════════════════════════════════════════════════════════════════════

  -- We use hasEigenvector_of_isMaxOn/isMinOn from Rayleigh.lean
  -- The key is to show the Rayleigh quotient attains its extremum for compact T

  -- Step 1: Define the Rayleigh quotient supremum and infimum
  let M_sup := ⨆ x : { x : E // x ≠ 0 }, T.rayleighQuotient x
  let M_inf := ⨅ x : { x : E // x ≠ 0 }, T.rayleighQuotient x

  -- Step 2: For self-adjoint T, ‖T‖ = max(|M_sup|, |M_inf|)
  -- So either |M_sup| = ‖T‖ or |M_inf| = ‖T‖

  -- Step 3: The closure of T applied to the unit sphere is compact
  have hK : IsCompact (closure (T '' Metric.sphere (0 : E) 1)) := by
    have h1 : Metric.sphere (0 : E) 1 ⊆ Metric.closedBall 0 1 := Metric.sphere_subset_closedBall
    have h2 : T '' Metric.sphere (0 : E) 1 ⊆ T '' Metric.closedBall 0 1 := Set.image_mono h1
    have h3 : closure (T '' Metric.sphere (0 : E) 1) ⊆ closure (T '' Metric.closedBall 0 1) :=
      closure_mono h2
    exact (hT.isCompact_closure_image_closedBall 1).of_isClosed_subset isClosed_closure h3

  -- Step 4: Get a maximizing sequence and use compactness
  -- Since T ≠ 0, the sphere is nonempty and there exist nonzero vectors
  have hE_nontrivial : Nontrivial E := by
    by_contra h
    have : Subsingleton E := not_nontrivial_iff_subsingleton.mp h
    have hT_eq : T = 0 := by ext x; simp [Subsingleton.eq_zero x]
    exact hT_ne hT_eq

  obtain ⟨v₀, hv₀⟩ : ∃ v₀ : E, v₀ ≠ 0 := exists_ne 0
  have hSphere_nonempty : (Metric.sphere (0 : E) 1).Nonempty := by
    use (‖v₀‖⁻¹ : 𝕜) • v₀
    simp only [Metric.mem_sphere, dist_zero_right, norm_smul]
    rw [norm_inv, RCLike.norm_coe_norm, inv_mul_cancel₀]
    exact norm_ne_zero_iff.mpr hv₀

  -- The proof proceeds by finding where the Rayleigh quotient achieves its maximum or minimum
  -- Using sequential compactness from the compact operator

  -- For the full proof, we need to:
  -- 1. Take a sequence xₙ in the unit sphere with ⟪Txₙ, xₙ⟫ approaching sup (or inf)
  -- 2. Extract a convergent subsequence of Txₙ using compactness
  -- 3. Show xₙ converges using the self-adjoint property
  -- 4. Apply hasEigenvector_of_isMaxOn or hasEigenvector_of_isMinOn

  -- The Rayleigh quotient is continuous
  have hRayleigh_cont : Continuous T.reApplyInnerSelf := T.reApplyInnerSelf_continuous

  -- For compact operators, we use sequential compactness
  -- The image of unit sphere under T is relatively compact

  -- Key technical step: Show that a maximizing sequence converges
  -- This requires showing ‖Txₙ - μxₙ‖ → 0 for the limiting eigenvalue μ

  -- We show existence by cases: either the sup or the inf gives us ‖T‖

  -- Case analysis: we consider whether sup ≥ |inf| or inf ≤ -|sup|
  -- In either case, we get an eigenvalue with absolute value ‖T‖

  -- The detailed proof requires careful analysis of the maximizing sequence.
  -- Here we provide the complete argument:

  -- First, note that for any unit vector x: |⟪Tx, x⟫| ≤ ‖T‖
  have hRayleigh_bound : ∀ x : E, ‖x‖ = 1 → |RCLike.re ⟪T x, x⟫| ≤ ‖T‖ := by
    intro x hx
    have h1 : ‖⟪T x, x⟫‖ ≤ ‖T x‖ * ‖x‖ := norm_inner_le_norm (T x) x
    have h2 : ‖T x‖ ≤ ‖T‖ * ‖x‖ := T.le_opNorm x
    calc |RCLike.re ⟪T x, x⟫| ≤ ‖⟪T x, x⟫‖ := RCLike.abs_re_le_norm _
      _ ≤ ‖T x‖ * ‖x‖ := h1
      _ ≤ ‖T‖ * ‖x‖ * ‖x‖ := by nlinarith [norm_nonneg x, norm_nonneg (T x)]
      _ = ‖T‖ := by rw [hx]; ring

  -- Local spectral-norm characterization for self-adjoint operators:
  -- the operator norm equals the supremum of the absolute Rayleigh quotient
  -- over the unit sphere.
  -- This is a standard result in functional analysis (see e.g., Reed-Simon Vol 1).
  have hOpNorm_eq_sSup_absRayleigh
      (hSA : IsSelfAdjoint T) :
      sSup (Set.range (fun x : {x : E // ‖x‖ = 1} ↦ |RCLike.re ⟪T x, x⟫|)) = ‖T‖ := by
    classical
    let S := Set.range (fun x : {x : E // ‖x‖ = 1} ↦ |RCLike.re ⟪T x, x⟫|)
    -- Upper bound: |⟪Tx, x⟫| ≤ ‖T‖ for unit vectors
    have hS_upper : ∀ r ∈ S, r ≤ ‖T‖ := by
      rintro r ⟨x, rfl⟩
      calc |RCLike.re ⟪T x, x⟫| ≤ ‖⟪T x, x⟫‖ := RCLike.abs_re_le_norm _
        _ ≤ ‖T (x : E)‖ * ‖(x : E)‖ := norm_inner_le_norm _ _
        _ ≤ ‖T‖ * ‖(x : E)‖ * ‖(x : E)‖ := by gcongr; exact T.le_opNorm _
        _ = ‖T‖ := by rw [x.property]; ring
    have hS_bdd : BddAbove S := ⟨‖T‖, hS_upper⟩
    have hS_nonneg : ∀ r ∈ S, 0 ≤ r := by
      rintro r ⟨x, rfl⟩
      exact abs_nonneg _
    -- Nonemptiness: we need a unit vector to show S is nonempty
    have hS_nonempty : S.Nonempty := by
      obtain ⟨v, hv⟩ := hSphere_nonempty
      have hv_norm : ‖v‖ = 1 := mem_sphere_zero_iff_norm.mp hv
      let u : {x : E // ‖x‖ = 1} := ⟨v, hv_norm⟩
      refine ⟨|RCLike.re ⟪T u, u⟫|, u, rfl⟩
    apply le_antisymm
    · exact csSup_le hS_nonempty hS_upper
    -- Lower bound: ‖T‖ ≤ sSup S using the Polarization Identity
    -- For self-adjoint T: 4 Re⟪Tx, y⟫ = Re⟪T(x+y), x+y⟫ - Re⟪T(x-y), x-y⟫
    -- Combined with the parallelogram law, this gives |Re⟪Tu, v⟫| ≤ M for unit u, v
    · let M := sSup S
      -- Step 1: M is an upper bound for |Re⟪Tu, u⟫| when ‖u‖ = 1
      have hM_bound : ∀ z : E, ‖z‖ = 1 → |RCLike.re ⟪T z, z⟫| ≤ M := fun z hz ↦
        le_csSup hS_bdd ⟨⟨z, hz⟩, rfl⟩
      -- Step 2: Extend to |Re⟪Tv, v⟫| ≤ M * ‖v‖² for all v
      have h_bound_all : ∀ v : E, |RCLike.re ⟪T v, v⟫| ≤ M * ‖v‖^2 := by
        intro v
        by_cases hv : v = 0
        · simp [hv]
        have hv_pos : 0 < ‖v‖ := norm_pos_iff.mpr hv
        let u := (↑‖v‖⁻¹ : 𝕜) • v
        have hu : ‖u‖ = 1 := by
          simp only [u, norm_smul, RCLike.norm_ofReal, abs_of_pos (inv_pos.mpr hv_pos)]
          exact inv_mul_cancel₀ (ne_of_gt hv_pos)
        have h_le := hM_bound u hu
        -- v = ‖v‖ • u, so ⟪Tv, v⟫ = ‖v‖² * ⟪Tu, u⟫
        have hv_eq : v = (↑‖v‖ : 𝕜) • u := by
          simp only [u, smul_smul, RCLike.ofReal_inv]
          rw [mul_inv_cancel₀ (RCLike.ofReal_ne_zero.mpr (ne_of_gt hv_pos)), one_smul]
        have h_inner_scale : ⟪T v, v⟫ = (↑(‖v‖^2) : 𝕜) * ⟪T u, u⟫ := by
          conv_lhs => rw [hv_eq]
          simp only [map_smul, inner_smul_left, inner_smul_right, RCLike.conj_ofReal]
          rw [sq, RCLike.ofReal_mul]; ring
        rw [h_inner_scale]
        simp only [RCLike.mul_re, RCLike.ofReal_re, RCLike.ofReal_im, zero_mul, sub_zero]
        rw [abs_mul, abs_of_pos (sq_pos_of_pos hv_pos)]
        calc ‖v‖^2 * |RCLike.re ⟪T u, u⟫| ≤ ‖v‖^2 * M :=
               mul_le_mul_of_nonneg_left h_le (sq_nonneg _)
          _ = M * ‖v‖^2 := mul_comm _ _
      -- Step 3: Self-adjoint symmetry: Re⟪Tx, y⟫ = Re⟪Ty, x⟫
      have hSA_symm : ∀ a b : E, RCLike.re ⟪T a, b⟫ = RCLike.re ⟪T b, a⟫ := by
        intro a b
        -- For self-adjoint T: ⟪Ta, b⟫ = ⟪a, Tb⟫ = conj ⟪Tb, a⟫
        -- Taking Re of both sides: Re⟪Ta, b⟫ = Re(conj ⟪Tb, a⟫) = Re⟪Tb, a⟫
        have h1 : ⟪T a, b⟫ = ⟪a, T b⟫ := by
          -- adjoint_inner_right: ⟪b, (adjoint T) a⟫ = ⟪T b, a⟫
          have h := ContinuousLinearMap.adjoint_inner_right T b a
          rw [hSA.adjoint_eq] at h
          -- h : ⟪b, T a⟫ = ⟪T b, a⟫
          -- We need: ⟪T a, b⟫ = ⟪a, T b⟫
          -- inner_conj_symm x y : conj ⟪y, x⟫ = ⟪x, y⟫
          -- So inner_conj_symm (T a) b : conj ⟪b, T a⟫ = ⟪T a, b⟫
          -- Thus (inner_conj_symm (T a) b).symm : ⟪T a, b⟫ = conj ⟪b, T a⟫
          have h2 : ⟪T a, b⟫ = starRingEnd 𝕜 ⟪b, T a⟫ := (inner_conj_symm (T a) b).symm
          have h3 : ⟪a, T b⟫ = starRingEnd 𝕜 ⟪T b, a⟫ := (inner_conj_symm a (T b)).symm
          rw [h2, h, h3]
        calc RCLike.re ⟪T a, b⟫ = RCLike.re ⟪a, T b⟫ := by rw [h1]
          _ = RCLike.re (starRingEnd 𝕜 ⟪T b, a⟫) := by rw [inner_conj_symm]
          _ = RCLike.re ⟪T b, a⟫ := RCLike.conj_re _
      -- Step 4: Polarization identity for self-adjoint T
      have h_polarization : ∀ x y : E,
          4 * RCLike.re ⟪T x, y⟫ = RCLike.re ⟪T (x+y), x+y⟫ - RCLike.re ⟪T (x-y), x-y⟫ := by
        intro x y
        simp only [map_add, map_sub, inner_add_left, inner_add_right,
          inner_sub_left, inner_sub_right, map_add, map_sub]
        have h_sym := hSA_symm x y
        linarith

[Levi G]

      -- Step 5: For unit vectors u, v: |Re⟪Tu, v⟫| ≤ M

      have h_unit_bound : ∀ u v : E, ‖u‖ = 1 → ‖v‖ = 1 → |RCLike.re ⟪T u, v⟫| ≤ M := by

        intro u v hu hv

        have h_pol := h_polarization u v

        -- Parallelogram law: ‖u+v‖² + ‖u-v‖² = 2(‖u‖² + ‖v‖²)

        have h_par : ‖u + v‖ * ‖u + v‖ + ‖u - v‖ * ‖u - v‖ = 2 * (‖u‖ * ‖u‖ + ‖v‖ * ‖v‖) :=

          parallelogram_law_with_norm 𝕜 u v

        rw [hu, hv] at h_par

        simp only [one_mul] at h_par

        -- Convert to ^2 form

        have h_par' : ‖u + v‖^2 + ‖u - v‖^2 = 4 := by

          simp only [sq]; linarith

        -- |4 * Re⟪Tu, v⟫| ≤ |Re⟪T(u+v), u+v⟫| + |Re⟪T(u-v), u-v⟫|

        have h_abs_le : |4 * RCLike.re ⟪T u, v⟫| ≤

            |RCLike.re ⟪T (u+v), u+v⟫| + |RCLike.re ⟪T (u-v), u-v⟫| := by

          rw [h_pol]

          exact abs_sub _ _

        calc |RCLike.re ⟪T u, v⟫|

          = |4 * RCLike.re ⟪T u, v⟫| / 4 := by

              rw [abs_mul, abs_of_pos (by norm_num : (0:ℝ) < 4)]; ring

          _ ≤ (|RCLike.re ⟪T (u+v), u+v⟫| + |RCLike.re ⟪T (u-v), u-v⟫|) / 4 := by

              apply div_le_div_of_nonneg_right h_abs_le (by norm_num : (0:ℝ) ≤ 4)

          _ ≤ (M * ‖u+v‖^2 + M * ‖u-v‖^2) / 4 := by

              apply div_le_div_of_nonneg_right _ (by norm_num : (0:ℝ) ≤ 4)

              exact add_le_add (h_bound_all _) (h_bound_all _)

          _ = M * (‖u+v‖^2 + ‖u-v‖^2) / 4 := by ring

          _ = M * 4 / 4 := by rw [h_par']

          _ = M := by ring

      -- Step 6: Use opNorm_le_bound with the polarization bound

      apply ContinuousLinearMap.opNorm_le_bound _ (Real.sSup_nonneg hS_nonneg)

      intro x

      by_cases hx : x = 0

      · simp [hx]

      have hx_pos : 0 < ‖x‖ := norm_pos_iff.mpr hx

      -- Normalize x to unit vector

      let u := (↑‖x‖⁻¹ : 𝕜) • x

      have hu : ‖u‖ = 1 := by

        simp only [u, norm_smul, RCLike.norm_ofReal, abs_of_pos (inv_pos.mpr hx_pos)]

        exact inv_mul_cancel₀ (ne_of_gt hx_pos)

      by_cases hTu : T u = 0

      · -- If Tu = 0, then Tx = ‖x‖ • Tu = 0, so ‖Tx‖ = 0 ≤ M * ‖x‖

        have hTx : T x = (↑‖x‖ : 𝕜) • T u := by

          simp only [u, map_smul, smul_smul, RCLike.ofReal_inv]

          rw [mul_inv_cancel₀ (RCLike.ofReal_ne_zero.mpr (ne_of_gt hx_pos)), one_smul]

        simp only [hTx, hTu, smul_zero, norm_zero]

        exact mul_nonneg (Real.sSup_nonneg hS_nonneg) (norm_nonneg _)

      -- Use v = Tu / ‖Tu‖ as the test vector

      have hTu_pos : 0 < ‖T u‖ := norm_pos_iff.mpr hTu

      let v := (↑‖T u‖⁻¹ : 𝕜) • T u

      have hv : ‖v‖ = 1 := by

        simp only [v, norm_smul, RCLike.norm_ofReal, abs_of_pos (inv_pos.mpr hTu_pos)]

        exact inv_mul_cancel₀ (ne_of_gt hTu_pos)

      -- ⟪Tu, v⟫ = ⟪Tu, Tu/‖Tu‖⟫ = ‖Tu‖

      have h_inner : RCLike.re ⟪T u, v⟫ = ‖T u‖ := by

        simp only [v, inner_smul_right]

        rw [inner_self_eq_norm_sq_to_K]

        simp only [RCLike.mul_re, RCLike.ofReal_re, RCLike.ofReal_im, zero_mul, sub_zero]

        have h_re : RCLike.re ((↑‖T u‖ : 𝕜) ^ 2) = ‖T u‖^2 := by

          rw [← RCLike.ofReal_pow]; exact RCLike.ofReal_re _

        rw [h_re]

        field_simp [ne_of_gt hTu_pos]

      have h_Tu_le : ‖T u‖ ≤ M := by

        have := h_unit_bound u v hu hv

        rw [h_inner, abs_of_pos hTu_pos] at this

        exact this

      -- Tx = ‖x‖ • Tu, so ‖Tx‖ = ‖x‖ * ‖Tu‖ ≤ ‖x‖ * M = M * ‖x‖

      have hTx : T x = (↑‖x‖ : 𝕜) • T u := by

        simp only [u, map_smul, smul_smul, RCLike.ofReal_inv]

        rw [mul_inv_cancel₀ (RCLike.ofReal_ne_zero.mpr (ne_of_gt hx_pos)), one_smul]

      calc ‖T x‖ = ‖(↑‖x‖ : 𝕜) • T u‖ := by rw [hTx]

        _ = ‖(↑‖x‖ : 𝕜)‖ * ‖T u‖ := norm_smul _ _

        _ = ‖x‖ * ‖T u‖ := by simp only [RCLike.norm_ofReal, abs_of_pos hx_pos]

        _ ≤ ‖x‖ * M := by nlinarith [h_Tu_le, hx_pos]

        _ = M * ‖x‖ := mul_comm _ _      -- Step 5: For unit vectors u, v: |Re⟪Tu, v⟫| ≤ M
      have h_unit_bound : ∀ u v : E, ‖u‖ = 1 → ‖v‖ = 1 → |RCLike.re ⟪T u, v⟫| ≤ M := by
        intro u v hu hv
        have h_pol := h_polarization u v
        -- Parallelogram law: ‖u+v‖² + ‖u-v‖² = 2(‖u‖² + ‖v‖²)
        have h_par : ‖u + v‖ * ‖u + v‖ + ‖u - v‖ * ‖u - v‖ = 2 * (‖u‖ * ‖u‖ + ‖v‖ * ‖v‖) :=
          parallelogram_law_with_norm 𝕜 u v
        rw [hu, hv] at h_par
        simp only [one_mul] at h_par
        -- Convert to ^2 form
        have h_par' : ‖u + v‖^2 + ‖u - v‖^2 = 4 := by
          simp only [sq]; linarith
        -- |4 * Re⟪Tu, v⟫| ≤ |Re⟪T(u+v), u+v⟫| + |Re⟪T(u-v), u-v⟫|
        have h_abs_le : |4 * RCLike.re ⟪T u, v⟫| ≤
            |RCLike.re ⟪T (u+v), u+v⟫| + |RCLike.re ⟪T (u-v), u-v⟫| := by
          rw [h_pol]
          exact abs_sub _ _
        calc |RCLike.re ⟪T u, v⟫|
          = |4 * RCLike.re ⟪T u, v⟫| / 4 := by
              rw [abs_mul, abs_of_pos (by norm_num : (0:ℝ) < 4)]; ring
          _ ≤ (|RCLike.re ⟪T (u+v), u+v⟫| + |RCLike.re ⟪T (u-v), u-v⟫|) / 4 := by
              apply div_le_div_of_nonneg_right h_abs_le (by norm_num : (0:ℝ) ≤ 4)
          _ ≤ (M * ‖u+v‖^2 + M * ‖u-v‖^2) / 4 := by
              apply div_le_div_of_nonneg_right _ (by norm_num : (0:ℝ) ≤ 4)
              exact add_le_add (h_bound_all _) (h_bound_all _)
          _ = M * (‖u+v‖^2 + ‖u-v‖^2) / 4 := by ring
          _ = M * 4 / 4 := by rw [h_par']
          _ = M := by ring
      -- Step 6: Use opNorm_le_bound with the polarization bound
      apply ContinuousLinearMap.opNorm_le_bound _ (Real.sSup_nonneg hS_nonneg)
      intro x
      by_cases hx : x = 0
      · simp [hx]
      have hx_pos : 0 < ‖x‖ := norm_pos_iff.mpr hx
      -- Normalize x to unit vector
      let u := (↑‖x‖⁻¹ : 𝕜) • x
      have hu : ‖u‖ = 1 := by
        simp only [u, norm_smul, RCLike.norm_ofReal, abs_of_pos (inv_pos.mpr hx_pos)]
        exact inv_mul_cancel₀ (ne_of_gt hx_pos)
      by_cases hTu : T u = 0
      · -- If Tu = 0, then Tx = ‖x‖ • Tu = 0, so ‖Tx‖ = 0 ≤ M * ‖x‖
        have hTx : T x = (↑‖x‖ : 𝕜) • T u := by
          simp only [u, map_smul, smul_smul, RCLike.ofReal_inv]
          rw [mul_inv_cancel₀ (RCLike.ofReal_ne_zero.mpr (ne_of_gt hx_pos)), one_smul]
        simp only [hTx, hTu, smul_zero, norm_zero]
        exact mul_nonneg (Real.sSup_nonneg hS_nonneg) (norm_nonneg _)
      -- Use v = Tu / ‖Tu‖ as the test vector
      have hTu_pos : 0 < ‖T u‖ := norm_pos_iff.mpr hTu
      let v := (↑‖T u‖⁻¹ : 𝕜) • T u
      have hv : ‖v‖ = 1 := by
        simp only [v, norm_smul, RCLike.norm_ofReal, abs_of_pos (inv_pos.mpr hTu_pos)]
        exact inv_mul_cancel₀ (ne_of_gt hTu_pos)
      -- ⟪Tu, v⟫ = ⟪Tu, Tu/‖Tu‖⟫ = ‖Tu‖
      have h_inner : RCLike.re ⟪T u, v⟫ = ‖T u‖ := by
        simp only [v, inner_smul_right]
        rw [inner_self_eq_norm_sq_to_K]
        simp only [RCLike.mul_re, RCLike.ofReal_re, RCLike.ofReal_im, zero_mul, sub_zero]
        have h_re : RCLike.re ((↑‖T u‖ : 𝕜) ^ 2) = ‖T u‖^2 := by
          rw [← RCLike.ofReal_pow]; exact RCLike.ofReal_re _
        rw [h_re]
        field_simp [ne_of_gt hTu_pos]
      have h_Tu_le : ‖T u‖ ≤ M := by
        have := h_unit_bound u v hu hv
        rw [h_inner, abs_of_pos hTu_pos] at this
        exact this
      -- Tx = ‖x‖ • Tu, so ‖Tx‖ = ‖x‖ * ‖Tu‖ ≤ ‖x‖ * M = M * ‖x‖
      have hTx : T x = (↑‖x‖ : 𝕜) • T u := by
        simp only [u, map_smul, smul_smul, RCLike.ofReal_inv]
        rw [mul_inv_cancel₀ (RCLike.ofReal_ne_zero.mpr (ne_of_gt hx_pos)), one_smul]
      calc ‖T x‖ = ‖(↑‖x‖ : 𝕜) • T u‖ := by rw [hTx]
        _ = ‖(↑‖x‖ : 𝕜)‖ * ‖T u‖ := norm_smul _ _
        _ = ‖x‖ * ‖T u‖ := by simp only [RCLike.norm_ofReal, abs_of_pos hx_pos]
        _ ≤ ‖x‖ * M := by nlinarith [h_Tu_le, hx_pos]
        _ = M * ‖x‖ := mul_comm _ _

[Levi G]

 -- Above we proved `hOpNorm_eq_sSup_absRayleigh`: the operator norm of a

  -- self-adjoint operator equals the supremum of the absolute Rayleigh quotient

  -- over the unit sphere. The proof uses the Polarization Identity:

  --   4 Re⟪Tx, y⟫ = Re⟪T(x+y), x+y⟫ - Re⟪T(x-y), x-y⟫

  -- Combined with the Parallelogram Law, this yields |Re⟪Tu, v⟫| ≤ M for

  -- unit vectors u, v (where M = sSup S), from which ‖T‖ ≤ M follows.

  --

  -- Now we use this to find eigenvectors. For compact T, we use the sequential

  -- characterization: there exists a sequence xₙ with ‖xₙ‖ = 1 and

  -- |⟪Txₙ, xₙ⟫| → ‖T‖. By compactness, Txₙ has a convergent subsequence, and

  -- the self-adjoint property gives convergence of xₙ itself.

  --

  -- These technical steps (maximizing sequence, subsequence extraction, and

  -- passage to the limit) are carried out explicitly below in

  -- `exists_eigenvector_from_maximizing_seq`, rather than being delegated

  -- directly to the Rayleigh lemmas `hasEigenvector_of_isMaxOn` and

  -- `hasEigenvector_of_isMinOn`.

  -- Helper lemma: from a maximizing sequence for the Rayleigh quotient,

  -- extract an eigenvector whose eigenvalue has absolute value equal to the

  -- limit of the Rayleigh values.

  have exists_eigenvector_from_maximizing_seq

      (x : ℕ → E) (hx_norm : ∀ n, ‖x n‖ = 1)

      (L : ℝ) (hL_nonneg : 0 ≤ L) (hL_eq_norm : L = ‖T‖)

      (hx_lim : Tendsto (fun n ↦ |RCLike.re ⟪T (x n), x n⟫|) atTop (𝓝 L)) :

      ∃ (μ : ℝ) (e : E), ‖e‖ = 1 ∧ T e = (μ : 𝕜) • e ∧ |μ| = L := by

    -- ═══════════════════════════════════════════════════════════════════════

    -- PROOF STRUCTURE:

    -- Step 1: Extract subsequence with constant sign using tendsto_abs_extract_subseq

    -- Step 2: Extract convergent subsequence of T x using compactness

    -- Step 3: Show x itself converges using tendsto_norm_sub_smul_zero

    -- Step 4: Pass to limit to get eigenvector equation

    -- ═══════════════════════════════════════════════════════════════════════

    -- Step 1: Extract subsequence where Re ⟪T x_n, x_n⟫ → σ * L (with |σ| = 1)

    obtain ⟨φ₁, σ, hφ₁_mono, hσ, hφ₁_lim⟩ :=

      tendsto_abs_extract_subseq hL_nonneg hx_lim

    -- The subsequence x ∘ φ₁ still has unit norm

    have hx_φ₁_norm : ∀ n, ‖x (φ₁ n)‖ = 1 := fun n ↦ hx_norm (φ₁ n)

    -- The eigenvalue candidate is μ = σ * L

    let μ := σ * L

    -- |μ| = |σ| * L = 1 * L = L (since |σ| = 1 and L ≥ 0)

    have hμ_abs : |μ| = L := by

      calc |μ| = |σ * L| := rfl

        _ = |σ| * |L| := abs_mul σ L

        _ = 1 * |L| := by rw [hσ]

        _ = |L| := one_mul _

        _ = L := abs_of_nonneg hL_nonneg

    -- Step 2: Use compactness to extract convergent subsequence of T (x (φ₁ n))

    -- The sequence T (x (φ₁ n)) lies in closure (T '' sphere 0 1) which is compact

    have hT_in_compact : ∀ n, T (x (φ₁ n)) ∈ closure (T '' Metric.sphere (0 : E) 1) := by

      intro n

      apply subset_closure

      exact ⟨x (φ₁ n), Metric.mem_sphere.mpr (by simp [hx_φ₁_norm n]), rfl⟩

    obtain ⟨y, hy_mem, φ₂, hφ₂_mono, hTx_lim⟩ := hK.tendsto_subseq hT_in_compact

    -- Combined subsequence

    let φ := φ₁ ∘ φ₂

    have hφ_mono : StrictMono φ := hφ₁_mono.comp hφ₂_mono

    -- Properties of the combined subsequence

    have hx_φ_norm : ∀ n, ‖x (φ n)‖ = 1 := fun n ↦ hx_norm (φ n)

    have hTx_φ_lim : Tendsto (fun n ↦ T (x (φ n))) atTop (𝓝 y) := hTx_lim

    have hRayleigh_φ_lim : Tendsto (fun n ↦ RCLike.re ⟪T (x (φ n)), x (φ n)⟫) atTop (𝓝 μ) :=

      hφ₁_lim.comp hφ₂_mono.tendsto_atTop

    -- Step 3: Show ‖T x_φ(n) - μ x_φ(n)‖ → 0 using tendsto_norm_sub_smul_zero

    have hμ_eq_norm : |μ| = ‖T‖ := by rw [hμ_abs, hL_eq_norm]

    have h_almost_eig : Tendsto (fun n ↦ ‖T (x (φ n)) - (μ : 𝕜) • x (φ n)‖) atTop (𝓝 0) :=

      IsSelfAdjoint.tendsto_norm_sub_smul_zero (T := T) hSA hx_φ_norm hμ_eq_norm hRayleigh_φ_lim

    -- Step 4: Since T x_φ(n) → y and ‖T x_φ(n) - μ x_φ(n)‖ → 0, we have μ x_φ(n) → y

    -- So x_φ(n) → y/μ (if μ ≠ 0)

    have h_smul_lim : Tendsto (fun n ↦ (μ : 𝕜) • x (φ n)) atTop (𝓝 y) := by

      -- T x_φ(n) - μ x_φ(n) → 0, so μ x_φ(n) → T x_φ(n) → y

      have h1 : Tendsto (fun n ↦ T (x (φ n)) - (μ : 𝕜) • x (φ n)) atTop (𝓝 0) :=

        tendsto_zero_iff_norm_tendsto_zero.mpr h_almost_eig

      have h2 := hTx_φ_lim.sub h1

      simp only [sub_zero] at h2

      -- h2 : Tendsto (fun n ↦ T (x (φ n)) - (T (x (φ n)) - μ • x (φ n))) atTop (𝓝 y)

      -- which simplifies to Tendsto (fun n ↦ μ • x (φ n)) atTop (𝓝 y)

      convert h2 using 1

      ext n

      simp only [sub_sub_cancel]

    -- Case split: L = 0 or L ≠ 0

    by_cases hL_zero : L = 0

    · -- If L = 0, then μ = σ * 0 = 0 and T = 0 (since ‖T‖ = L = 0)

      -- But T ≠ 0, contradiction

      exfalso

      have hT_zero : ‖T‖ = 0 := hL_eq_norm ▸ hL_zero

      rw [opNorm_zero_iff] at hT_zero

      exact hT_ne hT_zero

    · -- If L ≠ 0, then μ ≠ 0 and we can divide

      have hμ_ne : (μ : 𝕜) ≠ 0 := by

        simp only [ne_eq, RCLike.ofReal_eq_zero]

        intro hμ_zero

        have : |μ| = 0 := abs_eq_zero.mpr hμ_zero

        rw [hμ_abs] at this

        exact hL_zero this

      -- x_φ(n) = (1/μ) • (μ • x_φ(n)) → (1/μ) • y

      have hx_lim : Tendsto (fun n ↦ x (φ n)) atTop (𝓝 ((μ : 𝕜)⁻¹ • y)) := by

        have h := h_smul_lim.const_smul (μ : 𝕜)⁻¹

        simp only [smul_smul, inv_mul_cancel₀ hμ_ne, one_smul] at h

        exact h

      -- Let e = (1/μ) • y

      let e := (μ : 𝕜)⁻¹ • y

      -- e has unit norm by limit_of_unit_seq_is_unit

      have he_norm : ‖e‖ = 1 := limit_of_unit_seq_is_unit hx_φ_norm hx_lim

      -- T e = μ e: pass to the limit in T x_φ(n) = μ x_φ(n) + error

      have hTe : T e = (μ : 𝕜) • e := by

        -- T is continuous, so T x_φ(n) → T e

        have h_T_cont := T.continuous.tendsto e

        have hT_of_lim : Tendsto (fun n ↦ T (x (φ n))) atTop (𝓝 (T e)) :=

          h_T_cont.comp hx_lim

        -- But we also have T x_φ(n) → y, so T e = y

        have hTe_eq_y : T e = y := tendsto_nhds_unique hT_of_lim hTx_φ_lim

        -- And μ e = μ • (1/μ • y) = y

        have hμe_eq_y : (μ : 𝕜) • e = y := by

          calc (μ : 𝕜) • e = (μ : 𝕜) • ((μ : 𝕜)⁻¹ • y) := rfl

            _ = ((μ : 𝕜) * (μ : 𝕜)⁻¹) • y := by rw [smul_smul]

            _ = (1 : 𝕜) • y := by rw [mul_inv_cancel₀ hμ_ne]

            _ = y := one_smul _ _

        rw [hTe_eq_y, hμe_eq_y]

      exact ⟨μ, e, he_norm, hTe, hμ_abs⟩

  -- Lemma C: existence of an eigenvector with eigenvalue of maximal absolute value.

  -- We need to construct a maximizing sequence where |Re ⟪T x_n, x_n⟫| → ‖T‖

  have h_limit_exists : ∃ (μ : ℝ) (e : E), ‖e‖ = 1 ∧ T e = (μ : 𝕜) • e ∧ |μ| = ‖T‖ := by

    -- For self-adjoint operators, ‖T‖ = sup_{‖x‖=1} |Re ⟪T x, x⟫|

    -- We can construct a sequence approaching this supremum

    -- The set of Rayleigh quotient values is bounded above by ‖T‖

    let S : Set ℝ := Set.range (fun x : {x : E // ‖x‖ = 1} ↦ |RCLike.re ⟪T x, x⟫|)

    have hS_nonempty : S.Nonempty := by

      obtain ⟨v, hv⟩ := hSphere_nonempty

      have hv_norm : ‖v‖ = 1 := by

        rw [Metric.mem_sphere, dist_zero_right] at hv

        exact hv

      exact ⟨|RCLike.re ⟪T v, v⟫|, ⟨⟨v, hv_norm⟩, rfl⟩⟩

    have hS_bdd : BddAbove S := by

      use ‖T‖

      intro r ⟨x, hx⟩

      rw [← hx]

      exact hRayleigh_bound x x.property

    -- There exists a sequence approaching sSup S

    obtain ⟨u, hu_mono, hu_tendsto, hu_mem⟩ := exists_seq_tendsto_sSup hS_nonempty hS_bdd

    -- Extract the actual unit vectors from the sequence

    have h_unit_seq : ∃ (x : ℕ → E), (∀ n, ‖x n‖ = 1) ∧

        Tendsto (fun n ↦ |RCLike.re ⟪T (x n), x n⟫|) atTop (𝓝 (sSup S)) := by

      choose f hf using hu_mem

      use fun n ↦ (f n : E)

      constructor

      · intro n; exact (f n).property

      · convert hu_tendsto using 1

        ext n

        exact hf n

    obtain ⟨x, hx_norm, hx_lim⟩ := h_unit_seq

    -- The key fact: sSup S = ‖T‖ for self-adjoint operators

    -- This is the "Rayleigh quotient characterizes the norm" result

    have hsSup_eq_norm : sSup S = ‖T‖ := by

      -- This is exactly the content of the local spectral-norm lemma above.

      -- The definition of S matches the set used there.

      have h := hOpNorm_eq_sSup_absRayleigh hSA

      simpa [S] using h

    -- Now apply the helper lemma with L = ‖T‖

    have hL_nonneg : 0 ≤ ‖T‖ := norm_nonneg T

    rw [hsSup_eq_norm] at hx_lim

    exact exists_eigenvector_from_maximizing_seq x hx_norm ‖T‖ hL_nonneg rfl hx_lim

  exact h_limit_exists

[Aaron Liu]

This feels like way too long for a mathlib proof

[Levi G]

Aaron Liu said:

This feels like way too long for a mathlib proof

Thank you for your time spent looking at my proof for rayleigh quotient theorem for self-adjoint compact operators. Thank you for your concern. I get that it is too long for a mathlib proof however it plays a critical role in my spectral theorem proof. It does the heavy lifting in my spectral theorem where I use it as a lemma in proof by induction to get a base case then for the induction step in combination with the axiom of dependent choice to get a sequence of eigenvectors and eigenvalues.

[Yan Yablonovskiy 🇺🇦]

Levi G said:

Aaron Liu said:

This feels like way too long for a mathlib proof

Thank you for your time spent looking at my proof for rayleigh quotient theorem for self-adjoint compact operators. Thank you for your concern. I get that it is too long for a mathlib proof however it plays a critical role in my spectral theorem proof. It does the heavy lifting in my spectral theorem where I use it as a lemma in proof by induction to get a base case then for the induction step in combination with the axiom of dependent choice to get a sequence of eigenvectors and eigenvalues.

I think what everyone might like, is for you to at least try 'golf' and split up the proof and shorten it massively before sharing (which requires some manual effort and lean comprehension). You can look at mathlib conventions and current code adjacent to the theory.

It is a high cognitive load to review ungolfed and unidiomatic code, if you really want someone to help you and since its volunteering time, you gotta put in a bit more work to make it easier.

[Levi G]

@Yan Yablonovskiy 🇺🇦 Noted. Thank for you the feedback.

[Kevin Buzzard]

Alternatively, given that someone else is also working on this and they are writing idiomatic code because they know something about Lean, it might be best just to wait for their work (as has been pointed out already) rather than just continuing to generate unidiomatic AI code (which kind of feels like a waste of time to me, given the situation, as it will never be reviewed or merged).

[Niels Voss]

Sorry that my collaborators and I have not had much time to work on this recently; we are currently studying for finals and trying to complete another project by January. We are however willing to discuss design decisions in the meantime. The current status of the PR is that #32126 defines the singular values; however I will not personally consider this PR complete until we have proven that for finite-dimensional inner product spaces the singular values as defined by the approximation numbers are equal to the singular values as defined by the square roots of eigenvalues of T*T. We need to prove this before completing the PR because otherwise we risk the definition being incorrect.

We haven't had time to look too much into it, but we currently believe that the proper progression to this proof is the Von Neumann minimax theorem (#31560) implies the min-max theorem which implies the min-max principle for singular values which in turn implies the theorem we want.

Also, I think Kim Morrison's comment above makes sense:

I love the plan, but I hope that if there is not an implementation plan (i.e. a person saying they are doing it in a time frame) that this would not get in the way of someone contributing an implementation for matrices, which there is demand for and has been for a while. I think it should suffice in an implementation for matrices to include a module doc link to this thread, giving a pointer to the proposed upgrade path.

So I don't want to discourage people from working on the matrix SVD. However, I have very specific ideas for how the matrix SVD should be set up to maximize generality (see for instance my comment about merging the Full, Thin, Compact, and Truncated SVDs into a single definition). In particular, SVD has an above-average difficulty of even stating the theorem statement idiomatically in Lean, and there have been previous SVD formalizations which have been rejected simply due to unidiomatic definitions (e.g. #6042). Translating a proof that uses unidiomatic definitions to a proof that uses idiomatic definitions can actually be harder than starting from idiomatic definitions in the first place, which means that a proof with unidiomatic definitions might not count as any progress at all towards getting something into mathlib.

So I would request that anyone interested in working on the matrix SVD and getting it into mathlib discuss the theorem statement they intend to formalize along with all new definitions before attempting to work on the proof.

[Levi G]

Jireh Loreaux said:

Levi G the fact that it's a big file is part of the problem. At this point, you don't have experience contributing to Mathlib, so a contribution of this magnitude is very hard to review. If you start with smaller projects and learn good design patterns from receiving review, then it gets easier because your code is higher quality. Let me highlight one example as evidence:

theorem IsCompactOperator.eigenvalue_tendsto_zero
    (hT : IsCompactOperator T) (_hSA : IsSelfAdjoint T)
    {μ : ℕ → ℝ} {e : ℕ → E}
    (he_orth : Orthonormal 𝕜 e)
    (he_eigen : ∀ n, T (e n) = (μ n : 𝕜) • (e n))
    (hμ_antitone : Antitone (fun n ↦ |μ n|)) :
    Tendsto (fun n ↦ |μ n|) atTop (𝓝 0) := by

This theorem has nothing to do with selfadjiont operators, or real eigenvalues, or Antitone in any fundamental way. That is, what's really going on here is that since e is orthonormal, is converges weakly to zero, and compact operators take weakly convergent sequences to norm convergent ones.

Again, I would love for you to become a productive member of the community, but please start with smaller projects. It may even be the case that there are peices of your code worth salvaging, but I don't have the energy or the time to try to sift through 2000 lines of code looking for a diamond in the rough.

Does this address the issue in that example? It's part of my new code that I'll post later.

import Mathlib.Analysis.InnerProductSpace.Basic
import Mathlib.Analysis.Normed.Operator.Compact
import Mathlib.Topology.Algebra.InfiniteSum.Basic
import Mathlib.Analysis.InnerProductSpace.Adjoint

/-!
# Weak Convergence and Compact Operators

This file proves that eigenvalues of compact operators corresponding to orthonormal
sequences of eigenvectors tend to zero.

## Main results

* `Orthonormal.tendsto_inner_zero`: Orthonormal sequences converge weakly to zero.
* `IsCompactOperator.tendsto_of_tendsto_weakly_zero`: Compact operators map weakly
  convergent sequences to norm convergent sequences.
* `IsCompactOperator.eigenvalue_tendsto_zero`: The sequence of eigenvalues of a compact
  operator on a Hilbert space tends to zero (for orthonormal eigenvectors).
-/

open Filter Topology

variable {𝕜 : Type*} [RCLike 𝕜]
variable {E : Type*} [NormedAddCommGroup E] [InnerProductSpace 𝕜 E]

local notation "⟪" x ", " y "⟫" => @inner 𝕜 E _ x y

/-- **Orthonormal Sequences Converge Weakly to Zero** (Bessel's Inequality Consequence)

For any orthonormal sequence `(eₙ)` and any vector `x`, we have `⟪x, eₙ⟫ → 0` as `n → ∞`.
This follows from Bessel's inequality: `∑ |⟪x, eₙ⟫|² ≤ ‖x‖²`, which implies the terms
of the sum must tend to zero.
-/
theorem Orthonormal.tendsto_inner_zero {e : ℕ → E} (he : Orthonormal 𝕜 e) (x : E) :
    Tendsto (fun n ↦ ⟪x, e n⟫) atTop (𝓝 0) := by
  have h_summable : Summable (fun n ↦ ‖⟪e n, x⟫‖^2) := he.inner_products_summable x
  have h_sq_tendsto := h_summable.tendsto_atTop_zero
  have h_norm_tendsto : Tendsto (fun n ↦ ‖⟪e n, x⟫‖) atTop (𝓝 0) := by
    have h := h_sq_tendsto.sqrt
    simp only [Real.sqrt_sq (norm_nonneg _), Real.sqrt_zero] at h
    exact h
  have h_norm_swap : (fun n ↦ ‖⟪x, e n⟫‖) = (fun n ↦ ‖⟪e n, x⟫‖) := by
    funext n
    simp [norm_inner_symm]
  have h_norm_tendsto' : Tendsto (fun n ↦ ‖⟪x, e n⟫‖) atTop (𝓝 0) :=
    h_norm_swap ▸ h_norm_tendsto
  exact tendsto_zero_iff_norm_tendsto_zero.mpr h_norm_tendsto'

variable [CompleteSpace E]

/-- **Compact Operators are Weak-to-Norm Continuous on Bounded Sets**

If `T` is a compact operator and `(xₙ)` is a bounded sequence converging weakly to `0`,
then `T(xₙ) → 0` in norm.
-/
theorem IsCompactOperator.tendsto_of_tendsto_weakly_zero
    {T : E →L[𝕜] E} (hT : IsCompactOperator T) {x : ℕ → E}
    (hx_bdd : ∃ M : ℝ, ∀ n, ‖x n‖ ≤ M)
    (hx_weak : ∀ y : E, Tendsto (fun n ↦ ⟪y, x n⟫) atTop (𝓝 0)) :
    Tendsto (fun n ↦ T (x n)) atTop (𝓝 0) := by
  obtain ⟨M, hM⟩ := hx_bdd
  rcases hT.image_closedBall_subset_compact (f := (T : E →ₛₗ[RingHom.id 𝕜] E)) (M + 1) with
    ⟨K, hK_compact, hTK⟩
  have h_in_K : ∀ n, T (x n) ∈ K := fun n ↦ by
    apply hTK
    refine ⟨x n, ?_, rfl⟩
    simp only [Metric.mem_closedBall, dist_zero_right]
    calc ‖x n‖ ≤ M := hM n
      _ ≤ M + 1 := by linarith
  rw [Metric.tendsto_atTop]
  intro ε hε
  classical
  by_contra h_not_eventually
  have h_not_ev_lt : ¬ ∀ᶠ n in atTop, ‖T (x n)‖ < ε := by
    intro h_ev
    rcases (eventually_atTop.1 h_ev) with ⟨N, hN⟩
    exact h_not_eventually ⟨N, by intro n hn; simpa [dist_zero_right] using hN n hn⟩
  have h_freq : ∃ᶠ n in atTop, ε ≤ ‖T (x n)‖ :=
    (not_eventually.1 h_not_ev_lt).mono (by intro n hn; exact le_of_not_gt hn)
  obtain ⟨φ, hφ_mono, hφ_prop⟩ := extraction_of_frequently_atTop h_freq
  have h_in_K' : ∀ n, T (x (φ n)) ∈ K := fun n ↦ h_in_K (φ n)
  obtain ⟨y, _, ψ, hψ_mono, hTx_lim⟩ := hK_compact.tendsto_subseq h_in_K'
  have hy_zero : y = 0 := by
    have h_weak_y : ∀ z : E, ⟪z, y⟫ = 0 := fun z ↦ by
      have h1 : Tendsto (fun n ↦ ⟪z, T (x (φ (ψ n)))⟫) atTop (𝓝 ⟪z, y⟫) :=
        (Filter.Tendsto.inner tendsto_const_nhds hTx_lim)
      have h2 : Tendsto (fun n ↦ ⟪z, T (x (φ (ψ n)))⟫) atTop (𝓝 0) := by
        have h_adj : Tendsto (fun n ↦ ⟪(ContinuousLinearMap.adjoint T) z, x (φ (ψ n))⟫) atTop (𝓝 0) :=
          (hx_weak ((ContinuousLinearMap.adjoint T) z)).comp ((hφ_mono.comp hψ_mono).tendsto_atTop)
        simpa [ContinuousLinearMap.adjoint_inner_left] using h_adj
      exact tendsto_nhds_unique h1 h2
    exact inner_self_eq_zero.mp (h_weak_y y)
  have h_norm_lim : Tendsto (fun n ↦ ‖T (x (φ (ψ n)))‖) atTop (𝓝 0) := by
    have hTx_zero : Tendsto (fun n ↦ T (x (φ (ψ n)))) atTop (𝓝 0) := by simpa [hy_zero] using hTx_lim
    simpa [Function.comp, norm_zero] using (continuous_norm.tendsto (0 : E)).comp hTx_zero
  have h_ge_ε : ∀ n, ε ≤ ‖T (x (φ (ψ n)))‖ := fun n ↦ hφ_prop (ψ n)
  have h_eventually_small : ∀ᶠ n in atTop, ‖T (x (φ (ψ n)))‖ < ε := by
    rw [Metric.tendsto_atTop] at h_norm_lim
    obtain ⟨N, hN⟩ := h_norm_lim ε hε
    filter_upwards [Filter.Ici_mem_atTop N] with n hn
    simpa [dist_zero_right] using hN n hn
  obtain ⟨n, hn⟩ := h_eventually_small.exists
  exact (lt_of_lt_of_le hn (h_ge_ε n)).false

/-- Eigenvalues of a compact operator decay to zero along an orthonormal eigenbasis. -/
theorem IsCompactOperator.eigenvalue_tendsto_zero
    {T : E →L[𝕜] E} (hT : IsCompactOperator T)
    {μ : ℕ → 𝕜} {e : ℕ → E}
    (he_orth : Orthonormal 𝕜 e)
    (he_eigen : ∀ n, T (e n) = μ n • e n) :
    Tendsto (fun n ↦ ‖μ n‖) atTop (𝓝 0) := by
  have h_weak : ∀ y : E, Tendsto (fun n ↦ ⟪y, e n⟫) atTop (𝓝 0) :=
    fun y ↦ he_orth.tendsto_inner_zero y
  have h_bdd : ∃ M : ℝ, ∀ n, ‖e n‖ ≤ M := ⟨1, fun n ↦ le_of_eq (he_orth.norm_eq_one n)⟩
  have h_Te_zero : Tendsto (fun n ↦ T (e n)) atTop (𝓝 0) :=
    hT.tendsto_of_tendsto_weakly_zero h_bdd h_weak
  have h_norm_eq : ∀ n, ‖μ n‖ = ‖T (e n)‖ := fun n ↦ by
    rw [he_eigen n, norm_smul, he_orth.norm_eq_one, mul_one]
  simpa [h_norm_eq] using (continuous_norm.tendsto (0 : E)).comp h_Te_zero

[Michael Rothgang]

@Levi G In the interest of using my time well, I will not be reviewing your code on this matter, for all the reasons already mentioned up-thread. (Other maintainers may or may not do the same.) Make of that what you will, but I think it's fair for you to know.