-- Lean 4 proof step from Mathlib.Probability.Distributions.Gaussian.Multivariate (goal)
Goal:
ι : Type u_2
inst : DecidableEq ι
I J : Finset ι
μ : EuclideanSpace ℝ ↥I
S : Matrix ↥I ↥I ℝ
hS : S.PosSemidef
hJI : J ⊆ I
i j : ↥J
⊢ Eq (α := ℝ)
    ((((ContinuousLinearMap.toBilinForm
                  (covarianceBilin
                    (Measure.map (⇑(EuclideanSpace.restrict₂ hJI) : (a : EuclideanSpace ℝ ↥I) → EuclideanSpace ℝ ↥J)
                      (multivariateGaussian μ S))) :
                (a : EuclideanSpace ℝ ↥J) → EuclideanSpace ℝ ↥J →ₗ[ℝ] ℝ)
              ((OrthonormalBasis.toBasis (EuclideanSpace.basisFun ↥J ℝ) : (a : ↥J) → EuclideanSpace ℝ ↥J) i :
                EuclideanSpace ℝ ↥J) :
            EuclideanSpace ℝ ↥J →ₗ[ℝ] ℝ) :
          (a : EuclideanSpace ℝ ↥J) → ℝ)
        ((OrthonormalBasis.toBasis (EuclideanSpace.basisFun ↥J ℝ) : (a : ↥J) → EuclideanSpace ℝ ↥J) j :
          EuclideanSpace ℝ ↥J) :
      ℝ)
    ((((ContinuousLinearMap.toBilinForm
                  (covarianceBilin
                    (multivariateGaussian
                      ((EuclideanSpace.restrict₂ hJI : (a : EuclideanSpace ℝ ↥I) → EuclideanSpace ℝ ↥J) μ)
                      (S.submatrix (fun (i : ↥J) ↦ Subtype.mk (p := fun (x : ι) ↦ x ∈ I) (↑i : ι) ⋯) fun (i : ↥J) ↦
                        Subtype.mk (p := fun (x : ι) ↦ x ∈ I) (↑i : ι) ⋯))) :
                (a : EuclideanSpace ℝ ↥J) → EuclideanSpace ℝ ↥J →ₗ[ℝ] ℝ)
              ((OrthonormalBasis.toBasis (EuclideanSpace.basisFun ↥J ℝ) : (a : ↥J) → EuclideanSpace ℝ ↥J) i :
                EuclideanSpace ℝ ↥J) :
            EuclideanSpace ℝ ↥J →ₗ[ℝ] ℝ) :
          (a : EuclideanSpace ℝ ↥J) → ℝ)
        ((OrthonormalBasis.toBasis (EuclideanSpace.basisFun ↥J ℝ) : (a : ↥J) → EuclideanSpace ℝ ↥J) j :
          EuclideanSpace ℝ ↥J) :
      ℝ)
have hv : cov[(fun (u : EuclideanSpace ℝ ↥J) ↦
        ⟪(OrthonormalBasis.toBasis (EuclideanSpace.basisFun ↥J ℝ) : (a : ↥J) → EuclideanSpace ℝ ↥J) i, u⟫) ∘
      (⇑(EuclideanSpace.restrict₂ hJI) : (a : EuclideanSpace ℝ ↥I) → EuclideanSpace ℝ ↥J),
    (fun (u : EuclideanSpace ℝ ↥J) ↦
        ⟪(OrthonormalBasis.toBasis (EuclideanSpace.basisFun ↥J ℝ) : (a : ↥J) → EuclideanSpace ℝ ↥J) j, u⟫) ∘
      (⇑(EuclideanSpace.restrict₂ hJI) : (a : EuclideanSpace ℝ ↥I) → EuclideanSpace ℝ ↥J);
    multivariateGaussian μ S] =
  (((covarianceBilin
              (multivariateGaussian ((EuclideanSpace.restrict₂ hJI : (a : EuclideanSpace ℝ ↥I) → EuclideanSpace ℝ ↥J) μ)
                (S.submatrix (fun (i : ↥J) ↦ Subtype.mk (p := fun (x : ι) ↦ x ∈ I) (↑i : ι) sorry) fun (i : ↥J) ↦
                  Subtype.mk (p := fun (x : ι) ↦ x ∈ I) (↑i : ι) sorry)) :
            (a : EuclideanSpace ℝ ↥J) → EuclideanSpace ℝ ↥J →L[ℝ] ℝ)
          ((OrthonormalBasis.toBasis (EuclideanSpace.basisFun ↥J ℝ) : (a : ↥J) → EuclideanSpace ℝ ↥J) i :
            EuclideanSpace ℝ ↥J) :
        EuclideanSpace ℝ ↥J →L[ℝ] ℝ) :
      (a : EuclideanSpace ℝ ↥J) → ℝ)
    ((OrthonormalBasis.toBasis (EuclideanSpace.basisFun ↥J ℝ) : (a : ↥J) → EuclideanSpace ℝ ↥J) j : EuclideanSpace ℝ ↥J)
Tactic: rw [ContinuousLinearMap.toBilinForm_apply, ContinuousLinearMap.toBilinForm_apply,
      covarianceBilin_apply_eq_cov, covariance_map]