# Variation of Least Squares with Symmetric Positive Semi Definite (PSD) Constraint

I am trying to solve the following convex optimization problem: \begin{align} & \min_{W} && \sum_{i=1}^n (\mathbf{x}_{i}^TW\mathbf{x}_{i} - y_i)^2 \\\\ & s.t. && W \succcurlyeq 0 \\\\ & && W = W^T \end{align}

where $$\mathbf{x}_i \in \mathbb{R^p}$$, $$W \in \mathbb{R}^{p \times p}$$ and $$y_i \geq 0$$.

Without the positive semidefinite constraint, the problem is pretty straightforward. Requiring positive semidefiniteness, however, makes it a bit tricky.

I thought about using the fact that $$W \succcurlyeq 0$$ if and only if there exists a symmetric $$A$$ such that $$W = AA^T$$, and solving the equivalent problem

\begin{align} & \min_{A} && \sum_{i=1}^n (\mathbf{x}_{i}^TAA^T\mathbf{x}_{i} - y_i)^2 \\\\ &s.t. && A = A^T \end{align}

Letting $$a_{ij}$$ be the $$(i,j)th$$ element of A, this optimization function is quartic (fourth-order) with respect to the $$a_{ij}$$'s. Because of this, I am unsure of how to proceed.

I would be grateful if someone could point me in the right direction as to how to solve this problem.

• Haven't worked it through, but you might find it helpful to note that, if $W = A^\top A$, then$$x^\top_i W x_i = x^\top_i A^\top A x_i = (Ax_i)^\top (Ax_i) = \|Ax_i\|^2.$$ Apr 11 '20 at 0:01

When dealing with such form of matrix multiplications always remember the Vectorization Trick with Kronecker Product for Matrix Equations:

$${x}_{i}^{T} W {x}_{i} - {y}_{i} \Rightarrow \left({x}_{i}^{T} \otimes {x}_{i}^{T} \right) \operatorname{Vec} \left( W \right) - \operatorname{Vec} \left( {y}_{i} \right) = \left({x}_{i}^{T} \otimes {x}_{i}^{T} \right) \operatorname{Vec} \left( W \right) - {y}_{i}$$

Since the problem is given by summing over $${x}_{i}$$ one could build the matrix:

$$X = \begin{bmatrix} {x}_{1}^{T} \otimes {x}_{1}^{T} \\ {x}_{2}^{T} \otimes {x}_{2}^{T} \\ \vdots \\ {x}_{n}^{T} \otimes {x}_{n}^{T} \end{bmatrix}$$

Then:

$$\arg \min_{W} \sum_{i = 1}^{n} {\left( {x}_{i}^{T} W {x}_{i} - {y}_{i} \right)}^{2} = \arg \min_{W} {\left\| X \operatorname{Vec} \left( W \right) - \boldsymbol{y} \right\|}_{2}^{2}$$

Where $$\boldsymbol{y}$$ is the column vector composed by $${y}_{i}$$.

Now the above has nice form of regular Least Squares. The handling of the constraint can be done using Projected Gradient Descent Method. The projection onto the set of Symmetric Matrices and Positive Semi Definite (PSD) Matrices cone are given by:

1. $$\operatorname{Proj}_{\mathcal{S}^{n}} \left( A \right) = \frac{1}{2} \left( A + {A}^{T} \right)$$. See Orthogonal Projection of a Matrix onto the Set of Symmetric Matrices.
2. $$\operatorname{Proj}_{\mathcal{S}_{+}^{n}} \left( A \right) = Q {\Lambda}_{+} {Q}^{T}$$ where $$A = Q \Lambda {Q}^{T}$$ is the eigen decomposition of $$A$$ and $${\Lambda}_{+}$$ means we zero any negative values in $$\Lambda$$. See Find the Matrix Projection of a Symmetric Matrix onto the set of Symmetric Positive Semi Definite (PSD) Matrices.

Since both the Symmetric Matrices Set and PSD Cone are Linear Sub Space hence even the greedy iterative projection on the set will yield an orthogonal projection on the intersection of the 2 sets. See Orthogonal Projection onto the Intersection of Convex Sets.

So, with all the tools above one could create his own solver using basic tools with no need for external libraries (Which might be slow or not scale).

I implemented the Projected Gradient Descent Method with the above projections in MATLAB. I compared results to CVX to validate the solution. This is the solution: My implementation is vanilla Gradient Descent with constant Step Size and no acceleration. If you add those you'll see convergence which is order of magnitude faster (I guess few tens of iterations). Not bad for hand made solver.

The MATLAB Code is accessible in my StackExchange Mathematics Q3619669 GitHub Repository.

• This was a very clear and practical response. Thanks! Assuming I understand correctly, I have one question: what is the difference between iteratively projecting vs. projecting just once after finding the optimal solution to the unconstrained problem? Would they give different solutions? Apr 13 '20 at 15:43
• @tygaking, Please mark it as answer and +1. Regarding your question, you can't do that. Think of "Gradient Descent" as some kind of projection by itself. Hence what we do above is basically iterative projection onto iterative set, one for the gradient descent, one for each other.
– Royi
Apr 13 '20 at 15:45

Are you doing FGLS or something?

You could try substituting the constraint into the object. For the two by two case, for example, solve $$\sum_{i = 1}^N \left(x_i'\left[\array{ \array{w_{11} & w_{12}} \\ \array{w_{12}& w_{22} }} \right]x_i - y_i \right)^2$$ where $$w_{12} = w_{21}$$ now by construction. Then the matrix will be symmetric.

To ensure positive semi-definiteness, you can then use the standard principal minors test: $$w_{ii} \ge 0$$ for each $$i$$, $$w_{11} w_{22} - w_{12}^2 \ge 0$$, and so on, with the determinant of the upper left-hand minor weakly positive.

That at least subsumes positive semi-definiteness into concrete constraints and adjustments to the problem. That sounds like a nightmare to solve however, using Kuhn-Tucker. A simpler sufficient condition for semi-definiteness is the dominant diagonal condition, that $$w_{ii} \ge \sum_{j \neq i} |w_{ij}|$$ for each row $$i$$, which would be much more computationally tractable. Perhaps it could give you a good initial guess before you try relaxing it to the standard principal minors constraints.

• There is no need for this. While your approach might lead to problem formulated as Least Squares with Linear Equality and Inequality Constraints, it will also require some kind of iterative solver (As there is no closed form solution for LS with Inequality Constraints). Since we have the projection onto the Symmetric and PSD Matrices we can use it. See my answer. Still, nice idea of yours! Especially with the Diagonally Dominant Matrix. I might post a question and answer about this.
– Royi
Apr 12 '20 at 0:00

This is a convex optimization problem which can be easily formulated, and then numerically solved via a convex optimization tool, such as CVX, YALMIP, CVXPY, CVXR. This is a linear Semidefinite Programming problem (SDP), for which numerical solvers exist.

Here is the code for CVX.

Assume $$x_i$$ is the ith column of matrix $$X$$

cvx_begin
variable W(p,p) semidefinite % constrains W to be symmetric positive semidefinite
Objective = 0;
for i=1:n
Objective = Objective + square(X(:,i)'*W*X(:,i) - y(i))
end
minimize(Objective)
cvx_end


CVX will transform the problem into the form required by the solver, call the solver, and return the solution.

• I think the use of square_pos() here is wrong. It will lead to the wrong solution. You can use square() or if you want to be sure square_abs() or something like that.
– Royi
Apr 11 '20 at 23:35
• @Royi Right you are. Now fixed. Thanks for pointing it out. Do you want to take over for me as moderator at CVX Forum? Apr 11 '20 at 23:43
• Care to explain what you meant? Thanks...
– Royi
Apr 18 '20 at 22:43
• @Royi You can take over answering many of the questions there, and I can go into semi-retirement from that forum. Apr 18 '20 at 23:02
• I would answer you with a PM on CVX Forum. Yet the PM feature is disabled. Want to move it to an email?
– Royi
Apr 22 '20 at 15:08

This does seem to be a straight forward 'nonnegative' least squares though with the sdp constraint. n = 10; p = 5; X = zeros(n,p^2); for ii = 1:n x = randn(p,1); temp = xx'; X(ii,:) = temp(:)'; end y = randn(n,1); cvx_begin sdp variable W(p,p) semidefinite minimize(norm(XW(:)-y)) cvx_end