code: def linsys(x0, rho, P, q):
      
      return np.linalg.solve(rho * np.eye(q.shape[0]) + P, rho * x0.copy() + q)
docs: Proximal operator for the linear approximation Ax = b
  
      Minimizes the function:
  
      .. math:: f(x) = (1/2)||Ax-b||_2^2 = (1/2)x^TA^TAx - (b^TA)x + b^Tb
  
      Parameters
      ----------
      x0 : array_like
          The starting or initial point used in the proximal update step
  
      rho : float
          Momentum parameter for the proximal step (larger value -> stays closer to x0)
  
      P : array_like
          The symmetric matrix A^TA, where we are trying to approximate Ax=b
  
      q : array_like
          The vector A^Tb, where we are trying to approximate Ax=b
  
      Returns
      -------
      theta : array_like
          The parameter vector found after running the proximal update step