code: def cost_min2(self, alpha):
          
          n = self.V.dim()
          ax = alpha[:n]
          ay = alpha[n:]
  
          # ml = pyamg.ruge_stuben_solver(self.L)
          # # ml = pyamg.smoothed_aggregation_solver(self.L)
          # print(ml)
          # print()
          # print(self.L)
          # print()
          # x = ml.solve(ax, tol=1e-10)
          # print('residual: {}'.format(numpy.linalg.norm(ax - self.L*x)))
          # print()
          # print(ax)
          # print()
          # print(x)
          # exit(1)
  
          # x = sparse.linalg.spsolve(self.L, ax)
          # print('residual: {}'.format(numpy.linalg.norm(ax - self.L*x)))
          # exit(1)
  
          q2, r2 = self.get_q2_r2(ax, ay)
  
          Lax = self.L * ax
          Lay = self.L * ay
  
          out = [
              0.5 * numpy.dot(Lax, Lax),
              0.5 * numpy.dot(Lay, Lay),
              0.5 * numpy.dot(q2 - 1, q2 - 1),
              0.5 * numpy.dot(r2, r2),
          ]
  
          if self.num_f_eval % 10000 == 0:
              print("{:7d}     {:e} {:e} {:e} {:e}".format(self.num_f_eval, *out))
  
          self.num_f_eval += 1
          return numpy.sum(out)
docs: Residual formulation, Hessian is a low-rank update of the identity.