code: def gradient(self):
          
          L = self.L
          n = self.L.shape[0]
          grad = {"Lu": zeros((n, n, n * self._L.shape[1]))}
          for ii in range(self._L.shape[0] * self._L.shape[1]):
              row = ii // self._L.shape[1]
              col = ii % self._L.shape[1]
              grad["Lu"][row, :, ii] = L[:, col]
              grad["Lu"][:, row, ii] += L[:, col]
  
          return grad
docs: Derivative of the covariance matrix over the lower triangular, flat part of L.
  
          It is equal to
  
              ∂K/∂Lᵢⱼ = ALᵀ + LAᵀ,
  
          where Aᵢⱼ is an n×m matrix of zeros except at [Aᵢⱼ]ᵢⱼ=1.
  
          Returns
          -------
          Lu : ndarray
              Derivative of K over the lower-triangular, flat part of L.