glimix_core/cov/_kron2sum.py
def LdKL_dot(self, v, v1=None):
        """
        Implements L(∂K)Lᵀv.

        The array v can have one or two dimensions and the first dimension has to have
        size n⋅p.

        Let vec(V) = v. We have

            L(∂K)Lᵀ⋅v = ((Lₕ∂C₀Lₕᵀ) ⊗ (LₓGGᵀLₓᵀ))vec(V) = vec(LₓGGᵀLₓᵀVLₕ∂C₀Lₕᵀ),

        when the derivative is over the parameters of C₀. Similarly,

            L(∂K)Lᵀv = ((Lₕ∂C₁Lₕᵀ) ⊗ (LₓLₓᵀ))vec(V) = vec(LₓLₓᵀVLₕ∂C₁Lₕᵀ),

        over the parameters of C₁.
        """
        self._init_svd()

        def dot(a, b):
            r = tensordot(a, b, axes=([1], [0]))
            if a.ndim > b.ndim:
                return r.transpose([0, 2, 1])
            return r

        Lh = self.Lh
        V = unvec(v, (self.Lx.shape[0], -1) + v.shape[1:])
        LdKL_dot = {
            "C0.Lu": empty((v.shape[0],) + v.shape[1:] + (self._C0.Lu.shape[0],)),
            "C1.Lu": empty((v.shape[0],) + v.shape[1:] + (self._C1.Lu.shape[0],)),
        }

        dC0 = self._C0.gradient()["Lu"]
        for i in range(self._C0.Lu.shape[0]):
            t = dot(self._LxG, dot(self._LxG.T, dot(V, Lh @ dC0[..., i] @ Lh.T)))
            LdKL_dot["C0.Lu"][..., i] = t.reshape((-1,) + t.shape[2:], order="F")

        dC1 = self._C1.gradient()["Lu"]
        for i in range(self._C1.Lu.shape[0]):
            t = dot(V, Lh @ dC1[..., i] @ Lh.T)
            LdKL_dot["C1.Lu"][..., i] = t.reshape((-1,) + t.shape[2:], order="F")

        return LdKL_dot
def LdKL_dot(self, v, v1=None):
        """
        Implements L(∂K)Lᵀv.

        The array v can have one or two dimensions and the first dimension has to have
        size n⋅p.

        Let vec(V) = v. We have

            L(∂K)Lᵀ⋅v = ((Lₕ∂C₀Lₕᵀ) ⊗ (LₓGGᵀLₓᵀ))vec(V) = vec(LₓGGᵀLₓᵀVLₕ∂C₀Lₕᵀ),

        when the derivative is over the parameters of C₀. Similarly,

            L(∂K)Lᵀv = ((Lₕ∂C₁Lₕᵀ) ⊗ (LₓLₓᵀ))vec(V) = vec(LₓLₓᵀVLₕ∂C₁Lₕᵀ),

        over the parameters of C₁.
        """
        self._init_svd()

        def dot(a, b):
            r = tensordot(a, b, axes=([1], [0]))
            if a.ndim > b.ndim:
                return r.transpose([0, 2, 1])
            return r

        Lh = self.Lh
        V = unvec(v, (self.Lx.shape[0], -1) + v.shape[1:])
        LdKL_dot = {
            "C0.Lu": empty((v.shape[0],) + v.shape[1:] + (self._C0.Lu.shape[0],)),
            "C1.Lu": empty((v.shape[0],) + v.shape[1:] + (self._C1.Lu.shape[0],)),
        }

        dC0 = self._C0.gradient()["Lu"]
        for i in range(self._C0.Lu.shape[0]):
            t = dot(self._LxG, dot(self._LxG.T, dot(V, Lh @ dC0[..., i] @ Lh.T)))
            LdKL_dot["C0.Lu"][..., i] = t.reshape((-1,) + t.shape[2:], order="F")

        dC1 = self._C1.gradient()["Lu"]
        for i in range(self._C1.Lu.shape[0]):
            t = dot(V, Lh @ dC1[..., i] @ Lh.T)
            LdKL_dot["C1.Lu"][..., i] = t.reshape((-1,) + t.shape[2:], order="F")

        return LdKL_dot
Implements L(∂K)Lᵀv.

        The array v can have one or two dimensions and the first dimension has to have
        size n⋅p.

        Let vec(V) = v. We have

            L(∂K)Lᵀ⋅v = ((Lₕ∂C₀Lₕᵀ) ⊗ (LₓGGᵀLₓᵀ))vec(V) = vec(LₓGGᵀLₓᵀVLₕ∂C₀Lₕᵀ),

        when the derivative is over the parameters of C₀. Similarly,

            L(∂K)Lᵀv = ((Lₕ∂C₁Lₕᵀ) ⊗ (LₓLₓᵀ))vec(V) = vec(LₓLₓᵀVLₕ∂C₁Lₕᵀ),

        over the parameters of C₁.