title: cvxpy: best strategy for enforcing smoothness in second dimension for large problems
body: <p>There might be an answer for that floating around somewhere but I wasn't able to find it.</p>
  
  <p>I want to minimize <img src="https://chart.googleapis.com/chart?cht=tx&amp;chl=%5C%7CAX-B%5C%7C%5E2_F%20%2B%20%5Clambda%20%5C%7CXD%5ET%5C%7C%5E2_F" alt="formula"> with variable <code>X &gt;= 0</code> and  1st derivative matrix <code>D</code>, so that <code>X</code> is smooth in column-direction and with relatively large data. </p>
  
  <p>In the past I have used various ways of solving this (e.g. with <code>scipy.optimize.lsq_linear</code> or <a href="https://pylops.readthedocs.io/en/latest/" rel="nofollow noreferrer">pylops</a>) and now wanted to try out <a href="https://www.cvxpy.org/index.html" rel="nofollow noreferrer">cvxpy</a> with the generic approach below (using SCS, because the problem wouldn't fit into memory otherwise):</p>
  
  <pre><code>def cvxpy_solve_smooth(A, B, lambd=0):
      D = get_D(B.shape[1])
      X = cp.Variable((A.shape[1], B.shape[1]))
      cost1 = cp.sum_squares(A @ X - B)
      cost2 = lambd * cp.sum_squares(X @ D.T)
      constr = [0 &lt;= X]
      prob = cp.Problem(cp.Minimize(cost1 + cost2), constr)
      prob.solve(solver=cp.SCS)
      return X.value
  
  def get_D(n):
      D = np.diag(np.ones(n - 1), 1)
      np.fill_diagonal(D, -1)
      D = D[:-1]
      return D
  </code></pre>
  
  <p>However, this is considerably slower than <code>scipy.optimize.lsq_linear</code> with sparse matrices. What can I do in terms of problem formulation, solver options, cvxpy advanced features etc. to improve performance? </p>