If one wants to solve many linear systems with the same matrix, $$\mathbf A\mathbf x_1 = \mathbf b_1, \quad \mathbf A\mathbf x_2 = \mathbf b_2, \quad \ldots, \quad \mathbf A\mathbf x_k = \mathbf b_k,$$ it is preferable to first precompute a factorization of the matrix $\mathbf A$, so that each system with different right-hand sides can be solved efficiently. In my case, $\mathbf A$ is symmetric and positive definite, so the Cholesky factorization works well.

Of course, solving $\mathbf A\mathbf x_i = \mathbf b_i$ is equivalent to solving the unconstrained quadratic minimization $\min \tfrac12\mathbf x_i^T\mathbf A\mathbf x_i - \mathbf b_i^T\mathbf x_i$. I am interested in solving many such problems, but with additional inequality constraints, $$\begin{gather} \min & \tfrac12\mathbf x_i^T\mathbf A\mathbf x - \mathbf b_i^T\mathbf x_i \\ \text{s.t. } & \mathbf C_i\mathbf x \succeq \mathbf d_i. \end{gather}$$ The constraints and the linear term of the objective both vary between problems, but the quadratic term remains constant. Is there an analogous method for speeding up such problems by exploiting a factorization of $\mathbf A$? For equality constraints ($\mathbf C_i\mathbf x = \mathbf d_i$) we have tried the Uzawa method, but it does not support inequalities.

More details: The matrix $\mathbf A$ is sparse, symmetric, and positive definite. Its sparse Cholesky factorization is available, but I can also precompute other factorizations if necessary. If a general solution is not possible, I am also interested in the special case where each constraint acts on only one degree of freedom, i.e. each row of $\mathbf C_i$ has only one nonzero entry.


I don't give the details, only some ideas.

Let $f(x)=1/2x^TAx-b^Tx,g(x)=Cx-d=[g_1,\cdots,g_n]^T$ and $L_i$ be the rows of $C$.

Then the equations in the $2n$ real unknowns $x,\lambda_1,\cdots,\lambda_n$, are $Df(x)-\sum_i\lambda_i Dg_i(x)=0,\lambda_i g_i(x)=0$, that is to say $x^TA-b^T-\sum_i\lambda_i L_i=0,\lambda_i g_i(x)=0$. Finally, if $\Lambda=[\lambda_1,\cdots,\lambda_n]^T$, then the $2n$ linear equations are $$Ax=b+C^T\Lambda,\lambda_i g_i(x)=0.$$ Kuhn-Tucker gives more precise conditions...

Of course, there are $2^n$ choices for the column $\Lambda$.

For example, assume that $n=4$ and $\lambda_3=\lambda_4=0$. Then the matrix of the system in $x_1,x_2,x_3,x_4,\lambda_1,\lambda_2$ is directed by the $6\times 6$ symmetric matrix: $$\begin{pmatrix}A&L_1^T&L_2^T\\L_1&0&0\\L_2&0&0\end{pmatrix}.$$ Unfortunately, this matrix is not necessarily $>0$. Yet, we can proceed as follows:

$x=A^{-1}(-\lambda_1L_1^T-\lambda_2L_2^T+b)$ where $\lambda_1,\lambda_2$ satisfy the (corrected) system: $\lambda_1L_1A^{-1}L_1^T+\lambda_2L_1A^{-1}L_2^T-L_1A^{-1}b+d_1=0$


Finally, the key is to calculate $A^{-1}L_1^T,A^{-1}L_2^T,A^{-1}b$, that is easy, using the Choleski decomposition of $A$.

EDIT. I did not read the hypothesis: "$A$ symmetric $>0$"... Then $f$ is strictly convex and the set $S$ defined by $Cx\geq d$ is convex. Thus, if you obtain a local minimum on $S$, it is the unique global minimum on $S$. Moreover, if $x_0=A^{-1}b\notin S$, then the minimum is reached on the edge of $S$, that is, some $\lambda_i$'s are not $0$. Since the degree of $f$ is $2$, you can use QPSolve, a library in Maple or the equivalent in Matlab. Anyway, all available methods are iterative; in particular, choose for the initial point, a point on the edge of $S$ that is close to $x_0$.

  • $\begingroup$ You have formed the saddle point matrix $\begin{bmatrix}A & C^T \\ C & 0\end{bmatrix}$. However the second row yields the equation $Cx = d$, not $Cx \succeq d$. $\endgroup$ – Rahul Apr 20 '17 at 7:35
  • $\begingroup$ See my toy example where I corrected the last two equations. When $\lambda_j\not= 0$, the condition is $L_j x=d_j$; when $\lambda_j=0$, it is $L_j x>d_j$; at the end of the calculation, you check whether this last condition is satisfied (or you use KK conditions). Here $L_1 x>d_1,L_2 x>d_2,L_3 x=d_3,L_4 x=d_4$, that is, two equalities. $\endgroup$ – loup blanc Apr 20 '17 at 17:32
  • $\begingroup$ So you are proposing to solve for $\lambda$ assuming $Cx=d$, and then check if it satisfies the KKT conditions? There are two difficulties here: (i) Your last two equations form a system with matrix $CA^{-1}C^T$, the Schur complement of $A$ in the saddle point matrix. Forming that matrix is easy as you say, but solving for $\lambda$ is not. (ii) Even if you get $\lambda$ you have to check whether it satisfies the KKT conditions $\lambda_i\ge0$. If some components do not pass, do you remove the inactive constraints and repeat? Do you have a reference that proves convergence for this process? $\endgroup$ – Rahul Apr 20 '17 at 18:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.