docs: def updateD_H(self, x):
          """
          Compute Hessian for update of D
  
          See [2] for derivation of Hessian
          """
          self.precompute(x)
          H = zeros((len(x), len(x)))
          Ai = zeros(self.A.shape[0])
          Aj = zeros(Ai.shape)
          for i in range(len(x)):
              Ai = self.A[:, i]
              ti = dot(self.AD, outer(self.R[:, i], Ai)) + dot(outer(Ai, self.R[i, :]), self.ADt)
  
              for j in range(i, len(x)):
                  Aj = self.A[:, j]
                  tj = outer(Ai, Aj)
                  H[i, j] = (
                      self.E * (self.R[i, j] * tj + self.R[j, i] * tj.T) -
                      ti * (
                          dot(self.AD, outer(self.R[:, j], Aj)) +
                          dot(outer(Aj, self.R[j, :]), self.ADt)
                      )
                  ).sum()
                  H[j, i] = H[i, j]
          H *= -2
          e = eigvals(H).min()
          H = H + (eye(H.shape[0]) * e)
          return H