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.
 A: 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:


*

*$ \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.

*$ \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.
A: 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.
A: 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.
A: 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
