Problem formulation in CVX toolbox in Python

Question

I have an optimization problem to solve, which is a constrained least squares. I am trying to solve it using CVX toolbox in Python.

Here is the optimization problem:

$$\hat{H} = \arg \min_H \vert \vert A \: \text{vec}(H) - b\vert \vert_2^2$$ $$\text{subject to} \quad H^TH = I$$

Here $H \in \mathbb{R}^{p \times p}$ where $p = 2$.

Question is how to formulate it correctly in one of the CVX frameworks? Please bear in mind that I am comparatively new to CVX and still getting the hang of the syntax.

Given the quadratic constraint, I tried to do it this way

$$\hat{H} = \arg \min_H \: \textbf{tr}(I - H^TH)$$ $$\text{subject to} \quad A \: \text{vec}(H) = b$$

But this does not follow some DCP rules. The choice of using trace was to make sure the cost function is a scalar as required by CVX. How else can I formulate this? Please let me know if more clarifications are required!

Here is my code that I tried:

H1 = cp.Variable((nDim, nDim))
prob = cp.Problem(cp.Minimize(cp.trace(np.identity(nDim) - H1.T @ H1)),
                  [Phi @ cp.vec(H1) == cp.vec(B1)])
prob.solve()

If needed I can also update the code with some randomly generated matrices $A$ and $b$.

Answer 1

As noted in the comments , orthogonality is not a convex condition. Hence it cannot be solved by convex optimization and by CVX. Also finding exact solutions is unlikely.

However, here is a proposition of a probably non-optimal workaround partly using convex optimization. (To be tested)

First, search in the whole space. Solve the least square as usual

$$\min_H \|A~\text{vec}(H)-b\|_2^2$$

We note its solution $H_{LS}$, which is probably not orthonormal. Second, project this solution in the space of orthogonal matrices.

For this, perform its SVD

$$H_{LS}=USV^T$$

And finally, return $\hat H = UV^T$.

If you are lucky and the singular values of $S$ are both close to $1$, we can be pretty confident in the solution. Otherwise, you should probably make a sanity check somehow on the solution. Or use another method where the residual error is minimized "in the end" instead of minimizing then trying to search for a close solution in the set of orthogonal matrices (may be gradient descent with projection into orthogonal matrices using SVD as above at each step).

Update

If $H$ is only 2x2 like in your question, orthogonal matrices are rotations or reflexions. Both transformations are parametrized by $\theta \in[0,2\pi]$ in the plane so you can easily test a lot of $\theta$ for each type and keep the one minimizing at best the cost function (this also provides a sanity check for the above method). For larger dimensions it will quickly get out of hand though.