Just a suggestion. It seemed odd not to see this in the project. I’m not sure where best to insert it. The long-time contributors would know better.
import Mathlib
open Matrix
/-!
# The differential of the determinant
Jacobi's formula for square matrices of arbitrary finite size, including at
singular matrices. The algebraic linear functional is defined over a
commutative ring; the analytic results hold over any `RCLike` field, hence
uniformly for both `ℝ` and `ℂ`.
-/
namespace Matrix
variable {ι : Type*} [Fintype ι] [DecidableEq ι]
section Algebraic
variable {R : Type*} [CommRing R]
/-- The algebraic Jacobi functional at `M`:
`H ↦ ∑ i, j, adjugate M j i * H i j`.
The transposed indices agree with Mathlib's convention for `adjugate`: the
coefficient of the entry `(i, j)` is the cofactor `(i, j)`, namely the
`(j, i)` entry of the adjugate. -/
def detDifferentialLinearMap (M : Matrix ι ι R) :
Matrix ι ι R →ₗ[R] R where
toFun H := ∑ i : ι, ∑ j : ι, M.adjugate j i * H i j
map_add' H K := by
simp only [add_apply, mul_add, Finset.sum_add_distrib]
map_smul' c H := by
simp only [smul_eq_mul, RingHom.id_apply, smul_apply]
simp_rw [mul_left_comm (M.adjugate _ _) c]
simp only [← Finset.mul_sum]
@[simp] theorem detDifferentialLinearMap_apply (M H : Matrix ι ι R) :
detDifferentialLinearMap M H =
∑ i : ι, ∑ j : ι, M.adjugate j i * H i j := rfl
end Algebraic
section Analytic
variable {𝕜 : Type*} [RCLike 𝕜]
noncomputable local instance matrixNormedAddCommGroup :
NormedAddCommGroup (Matrix ι ι 𝕜) := Matrix.normedAddCommGroup
noncomputable local instance matrixNormedSpace :
NormedSpace 𝕜 (Matrix ι ι 𝕜) := Matrix.normedSpace
/-- The continuous Jacobi functional at `M`. Since the matrix space is
finite-dimensional, the algebraic Jacobi functional is automatically continuous. -/
noncomputable def detDifferential (M : Matrix ι ι 𝕜) :
Matrix ι ι 𝕜 →L[𝕜] 𝕜 :=
(detDifferentialLinearMap M).toContinuousLinearMap
@[simp] theorem detDifferential_apply (M H : Matrix ι ι 𝕜) :
detDifferential M H =
∑ i : ι, ∑ j : ι, M.adjugate j i * H i j := rfl
/-- Jacobi's formula without an invertibility assumption. -/
theorem hasStrictFDerivAt_det_adjugate (M : Matrix ι ι 𝕜) :
HasStrictFDerivAt (fun A : Matrix ι ι 𝕜 => A.det)
(detDifferential M) M := by
-- The determinant of a matrix is a polynomial function in its entries, hence it is infinitely differentiable.
have h_det_smooth : ContDiff 𝕜 ⊤ (fun A : (ι → ι → 𝕜) => Matrix.det (Matrix.of (fun i j => A i j))) := by
simp +decide [Matrix.det_apply']
fun_prop (disch := norm_num)
have h_det_deriv : ∀ (v : ι → ι → 𝕜), deriv (fun t : 𝕜 => Matrix.det (Matrix.of (fun i j => M i j + t * v i j))) 0 = ∑ i, ∑ j, Matrix.adjugate M j i * v i j := by
simp +decide [Matrix.det_apply', Matrix.adjugate_apply, Matrix.det_apply'];
intro v
have h_det_deriv : deriv (fun t : 𝕜 => ∑ x : Equiv.Perm ι, (Equiv.Perm.sign x : 𝕜) * ∏ i : ι, (M (x i) i + t * v (x i) i)) 0 = ∑ x : Equiv.Perm ι, (Equiv.Perm.sign x : 𝕜) * ∑ i : ι, (∏ j ∈ Finset.univ.erase i, M (x j) j) * v (x i) i := by
convert HasDerivAt.deriv _ using 1
have h_deriv : ∀ x : Equiv.Perm ι, HasDerivAt (fun t : 𝕜 => ∏ i : ι, (M (x i) i + t * v (x i) i)) (∑ i : ι, (∏ j ∈ Finset.univ.erase i, M (x j) j) * v (x i) i) 0 := by
intro x
have h_deriv : ∀ (f : ι → 𝕜 → 𝕜), (∀ i, HasDerivAt (f i) (v (x i) i) 0) → HasDerivAt (fun t => ∏ i, f i t) (∑ i, (∏ j ∈ Finset.univ.erase i, f j 0) * v (x i) i) 0 := by
exact fun f hf => HasDerivAt.fun_finsetProd fun i _ => hf i
simpa using h_deriv (fun i t => M (x i) i + t * v (x i) i) fun i => by simpa using HasDerivAt.const_add (M (x i) i) (HasDerivAt.mul (hasDerivAt_id 0) (hasDerivAt_const _ _))
convert HasDerivAt.sum fun x _ => HasDerivAt.const_mul (Equiv.Perm.sign x : 𝕜) (h_deriv x) using 1
exact
Eq.symm
(Finset.sum_fn Finset.univ fun c y =>
↑↑(Equiv.Perm.sign c) * ∏ i, (M (c i) i + y * v (c i) i))
rw [h_det_deriv, Finset.sum_comm]
simp +decide [Matrix.updateRow_apply, Finset.mul_sum _ _ _, mul_assoc, mul_comm]
refine' Finset.sum_comm.trans (Finset.sum_congr rfl fun i _ => _)
rw [Finset.sum_comm, Finset.sum_congr rfl]
intros
simp only [Finset.prod_ite]
rw [Finset.sum_eq_single (‹Equiv.Perm ι› i)] <;>
simp +decide [Finset.filter_ne', Finset.filter_eq']
intro j hj
simp only [Finset.prod_eq_zero_iff, Pi.single_apply]
exact Or.inr <| Or.inl <| ⟨ ‹Equiv.Perm ι›.symm j, by simp +decide, by simpa [Equiv.symm_apply_eq] using hj ⟩
refine' (h_det_smooth.contDiffAt.hasStrictFDerivAt _).congr_fderiv _
· decide +revert
· -- By definition of the derivative, we know that
have h_deriv : ∀ (v : ι → ι → 𝕜), deriv (fun t : 𝕜 => Matrix.det (Matrix.of (fun i j => M i j + t * v i j))) 0 = (fderiv 𝕜 (fun A : (ι → ι → 𝕜) => Matrix.det (Matrix.of (fun i j => A i j))) M) v := by
intro v
rw [deriv]
rw [show (fun t : 𝕜 => Matrix.det (Matrix.of (fun i j => M i j + t * v i j))) = (fun A : (ι → ι → 𝕜) => Matrix.det (Matrix.of (fun i j => A i j))) ∘ (fun t : 𝕜 => fun i j => M i j + t * v i j) by ext; rfl, fderiv_comp] <;> norm_num [h_det_smooth.contDiffAt.differentiableAt]
rw [deriv_pi] <;> norm_num [Finset.sum_add_distrib, mul_comm]
congr ; ext i j
rw [deriv_pi] <;> norm_num [Finset.sum_add_distrib, mul_comm]
exact ContinuousLinearMap.ext fun v => by
exact (h_deriv v).symm.trans ((h_det_deriv v).trans (detDifferential_apply M v).symm)
/-- The Fréchet derivative of the determinant as a continuous linear map. -/
theorem fderiv_det (M : Matrix ι ι 𝕜) :
fderiv 𝕜 (fun A : Matrix ι ι 𝕜 => A.det) M = detDifferential M := by
exact (hasStrictFDerivAt_det_adjugate M).hasFDerivAt.fderiv
/-- Jacobi's formula evaluated in a direction `H`. -/
theorem fderiv_det_apply (M H : Matrix ι ι 𝕜) :
fderiv 𝕜 (fun A : Matrix ι ι 𝕜 => A.det) M H =
∑ i : ι, ∑ j : ι, M.adjugate j i * H i j := by
rw [fderiv_det, detDifferential_apply]
end Analytic
end Matrix
Just a suggestion. It seemed odd not to see this in the project. I’m not sure where best to insert it. The long-time contributors would know better.