def _lml_gradient(self):
from scipy.linalg import cho_solve<EOL>terms = self._terms<EOL>dC0 = self._cov.C0.gradient()["<STR_LIT>"]<EOL>dC1 = self._cov.C1.gradient()["<STR_LIT>"]<EOL>b = terms["<STR_LIT:b>"]<EOL>W = terms["<STR_LIT>"]<EOL>Lh = terms["<STR_LIT>"]<EOL>Lz = terms["<STR_LIT>"]<EOL>WA = terms["<STR_LIT>"]<EOL>WL0 = terms["<STR_LIT>"]<EOL>YW = terms["<STR_LIT>"]<EOL>MRiM = terms["<STR_LIT>"]<EOL>MRiy = terms["<STR_LIT>"]<EOL>XRiM = terms["<STR_LIT>"]<EOL>XRiy = terms["<STR_LIT>"]<EOL>ZiXRiM = terms["<STR_LIT>"]<EOL>ZiXRiy = terms["<STR_LIT>"]<EOL>WdC0 = _mdot(W, dC0)<EOL>WdC1 = _mdot(W, dC1)<EOL>AWdC0 = _mdot(WA.T, dC0)<EOL>AWdC1 = _mdot(WA.T, dC1)<EOL>MR0M = _mkron(_mdot(AWdC0, WA), self._XGGX)<EOL>MR1M = _mkron(_mdot(AWdC1, WA), self._XX)<EOL>MR0X = _mkron(_mdot(AWdC0, WL0), self._XGGG)<EOL>MR1X = _mkron(_mdot(AWdC1, WL0), self._GX.T)<EOL>MR0y = vec(_mdot(self._XGGY, _mdot(WdC0, WA)))<EOL>MR1y = vec(_mdot(self._XY, WdC1, WA))<EOL>XR0X = _mkron(_mdot(WL0.T, dC0, WL0), self._GGGG)<EOL>XR1X = _mkron(_mdot(WL0.T, dC1, WL0), self._GG)<EOL>XR0y = vec(_mdot(self._GGGY, WdC0, WL0))<EOL>XR1y = vec(_mdot(self._GY, WdC1, WL0))<EOL>yR0y = vec(_mdot(self._GY, WdC0)).T @ vec(self._GY @ W)<EOL>yR1y = (YW.T * _mdot(self._Y, WdC1).T).T.sum(axis=(<NUM_LIT:0>, <NUM_LIT:1>))<EOL>ZiXR0X = cho_solve(Lz, XR0X)<EOL>ZiXR1X = cho_solve(Lz, XR1X)<EOL>ZiXR0y = cho_solve(Lz, XR0y)<EOL>ZiXR1y = cho_solve(Lz, XR1y)<EOL>MK0y = MR0y - _mdot(XRiM.T, ZiXR0y) - _mdot(MR0X, ZiXRiy)<EOL>MK0y += _mdot(XRiM.T, ZiXR0X, ZiXRiy)<EOL>MK1y = MR1y - _mdot(XRiM.T, ZiXR1y) - _mdot(MR1X, ZiXRiy)<EOL>MK1y += _mdot(XRiM.T, ZiXR1X, ZiXRiy)<EOL>yK0y = yR0y - <NUM_LIT:2> * XR0y.T @ ZiXRiy + ZiXRiy.T @ _mdot(XR0X, ZiXRiy)<EOL>yK1y = yR1y - <NUM_LIT:2> * XR1y.T @ ZiXRiy + ZiXRiy.T @ _mdot(XR1X, ZiXRiy)<EOL>MR0XZiXRiM = _mdot(MR0X, ZiXRiM)<EOL>MK0M = MR0M - MR0XZiXRiM - MR0XZiXRiM.transpose([<NUM_LIT:1>, <NUM_LIT:0>, <NUM_LIT:2>])<EOL>MK0M += _mdot(ZiXRiM.T, XR0X, ZiXRiM)<EOL>MR1XZiXRiM = _mdot(MR1X, ZiXRiM)<EOL>MK1M = MR1M - MR1XZiXRiM - MR1XZiXRiM.transpose([<NUM_LIT:1>, <NUM_LIT:0>, <NUM_LIT:2>])<EOL>MK1M += _mdot(ZiXRiM.T, XR1X, ZiXRiM)<EOL>MK0m = _mdot(MK0M, b)<EOL>mK0y = b.T @ MK0y<EOL>mK0m = b.T @ MK0m<EOL>MK1m = _mdot(MK1M, b)<EOL>mK1y = b.T @ MK1y<EOL>mK1m = b.T @ MK1m<EOL>XRim = XRiM @ b<EOL>MRim = MRiM @ b<EOL>db = {"<STR_LIT>": cho_solve(Lh, MK0m - MK0y), "<STR_LIT>": cho_solve(Lh, MK1m - MK1y)}<EOL>grad = {<EOL>"<STR_LIT>": -trace(WdC0) * self._trGG + trace(ZiXR0X),<EOL>"<STR_LIT>": -trace(WdC1) * self.nsamples + trace(ZiXR1X),<EOL>}<EOL>if self._restricted:<EOL><INDENT>grad["<STR_LIT>"] += cho_solve(Lh, MK0M).diagonal().sum(<NUM_LIT:1>)<EOL>grad["<STR_LIT>"] += cho_solve(Lh, MK1M).diagonal().sum(<NUM_LIT:1>)<EOL><DEDENT>mKiM = MRim.T - XRim.T @ ZiXRiM<EOL>yKiM = MRiy.T - XRiy.T @ ZiXRiM<EOL>grad["<STR_LIT>"] += yK0y - <NUM_LIT:2> * mK0y + mK0m - <NUM_LIT:2> * _mdot(mKiM, db["<STR_LIT>"])<EOL>grad["<STR_LIT>"] += <NUM_LIT:2> * _mdot(yKiM, db["<STR_LIT>"])<EOL>grad["<STR_LIT>"] += yK1y - <NUM_LIT:2> * mK1y + mK1m - <NUM_LIT:2> * _mdot(mKiM, db["<STR_LIT>"])<EOL>grad["<STR_LIT>"] += <NUM_LIT:2> * _mdot(yKiM, db["<STR_LIT>"])<EOL>grad["<STR_LIT>"] /= <NUM_LIT:2><EOL>grad["<STR_LIT>"] /= <NUM_LIT:2><EOL>return grad<EOL>
Gradient of the log of the marginal likelihood.

Let 𝐲 = vec(Y), 𝕂 = K⁻¹∂(K)K⁻¹, and H = MᵀK⁻¹M. The gradient is given by::

    2⋅∂log(p(𝐲)) = -tr(K⁻¹∂K) - tr(H⁻¹∂H) + 𝐲ᵀ𝕂𝐲 - 𝐦ᵀ𝕂(2⋅𝐲-𝐦)
        - 2⋅(𝐦-𝐲)ᵀK⁻¹∂(𝐦).

Observe that

    ∂𝛃 = -H⁻¹(∂H)𝛃 - H⁻¹Mᵀ𝕂𝐲 and ∂H = -Mᵀ𝕂M.

Let Z = I + XᵀR⁻¹X and 𝓡 = R⁻¹(∂K)R⁻¹. We use Woodbury matrix identity to
write ::

    𝐲ᵀ𝕂𝐲 = 𝐲ᵀ𝓡𝐲 - 2(𝐲ᵀ𝓡X)Z⁻¹(XᵀR⁻¹𝐲) + (𝐲ᵀR⁻¹X)Z⁻¹(Xᵀ𝓡X)Z⁻¹(XᵀR⁻¹𝐲)
    Mᵀ𝕂M = Mᵀ𝓡M - 2(Mᵀ𝓡X)Z⁻¹(XᵀR⁻¹M) + (MᵀR⁻¹X)Z⁻¹(Xᵀ𝓡X)Z⁻¹(XᵀR⁻¹M)
    Mᵀ𝕂𝐲 = Mᵀ𝓡𝐲 - (MᵀR⁻¹X)Z⁻¹(Xᵀ𝓡𝐲) - (Mᵀ𝓡X)Z⁻¹(XᵀR⁻¹𝐲)
          + (MᵀR⁻¹X)Z⁻¹(Xᵀ𝓡X)Z⁻¹(XᵀR⁻¹𝐲)
    H⁻¹   = MᵀR⁻¹M - (MᵀR⁻¹X)Z⁻¹(XᵀR⁻¹M),

where we have used parentheses to separate expressions
that we will compute separately. For example, we have ::

    𝐲ᵀ𝓡𝐲 = 𝐲ᵀ(U₁S₁⁻¹U₁ᵀ ⊗ I)(∂C₀ ⊗ GGᵀ)(U₁S₁⁻¹U₁ᵀ ⊗ I)𝐲
          = 𝐲ᵀ(U₁S₁⁻¹U₁ᵀ∂C₀ ⊗ G)(U₁S₁⁻¹U₁ᵀ ⊗ Gᵀ)𝐲
          = vec(GᵀYU₁S₁⁻¹U₁ᵀ∂C₀)ᵀvec(GᵀYU₁S₁⁻¹U₁ᵀ),

when the derivative is over the parameters of C₀. Otherwise, we have

    𝐲ᵀ𝓡𝐲 = vec(YU₁S₁⁻¹U₁ᵀ∂C₁)ᵀvec(YU₁S₁⁻¹U₁ᵀ).

The above equations can be more compactly written as

    𝐲ᵀ𝓡𝐲 = vec(EᵢᵀYW∂Cᵢ)ᵀvec(EᵢᵀYW),

where W = U₁S₁⁻¹U₁ᵀ, E₀ = G, and E₁ = I. We will now just state the results for
the other instances of the aBc form, which follow similar derivations::

    Xᵀ𝓡X = (L₀ᵀW∂CᵢWL₀) ⊗ (GᵀEᵢEᵢᵀG)
    Mᵀ𝓡y = (AᵀW∂Cᵢ⊗XᵀEᵢ)vec(EᵢᵀYW) = vec(XᵀEᵢEᵢᵀYW∂CᵢWA)
    Mᵀ𝓡X = AᵀW∂CᵢWL₀ ⊗ XᵀEᵢEᵢᵀG
    Mᵀ𝓡M = AᵀW∂CᵢWA ⊗ XᵀEᵢEᵢᵀX
    Xᵀ𝓡𝐲 = GᵀEᵢEᵢᵀYW∂CᵢWL₀

From Woodbury matrix identity and Kronecker product properties we have ::

    tr(K⁻¹∂K) = tr[W∂Cᵢ]tr[EᵢEᵢᵀ] - tr[Z⁻¹(Xᵀ𝓡X)]
    tr(H⁻¹∂H) = - tr[(MᵀR⁻¹M)(Mᵀ𝕂M)] + tr[(MᵀR⁻¹X)Z⁻¹(XᵀR⁻¹M)(Mᵀ𝕂M)]

Note also that ::

    ∂𝛃 = H⁻¹Mᵀ𝕂M𝛃 - H⁻¹Mᵀ𝕂𝐲.

Returns
-------
C0.Lu : ndarray
    Gradient of the log of the marginal likelihood over C₀ parameters.
C1.Lu : ndarray
    Gradient of the log of the marginal likelihood over C₁ parameters.