What method/algorithm for constrained multi-target regression

After some further research and inspiration I have figured out the following method:

I am using the vector notation of my problem with $underline{y}=begin{pmatrix}y_1y_2y_3end{pmatrix},~underline{x}=begin{pmatrix}1x_1x_2x_3end{pmatrix},~underline{underline{B}}=begin{pmatrix}beta_{10}&beta_{11}&beta_{12}&beta_{13} beta_{20}&beta_{21}&beta_{22}&beta_{23} beta_{30}&beta_{31}&beta_{32}&beta_{33} end{pmatrix}$ This results in the optimisation problem with $N$ examples

$underline{underline{hat{B}}}=underset{beta in C}{operatorname{argmin}}L(underline{underline{B}})=underset{beta in C}{operatorname{argmin}}sum_{i=1}^{N}lVert underline{y}_i - underline{underline{B}}cdotunderline{x}_i rVert^2,~~~~text{with}~~~~C= lbrace underline{underline{B}}in mathbb{R}^{3times4} mid beta_{11}+beta_{22}+beta_{33}=0rbrace$

For solving this, I am using projected gradient descent where $underline{underline{B}}$ for every iteration is projecet onto $C$. Normally the Projection operator would require another optimisation process $operatorname{proj}_C(x)=operatorname{argmin}_u biglbrace chi_C(u) + frac{1}{2}lVert u-xrVert^2bigrbrace$ with $chi_C$ the characteristic function.

Now comes the neat part: Since the constraint $C$ describes a flat subspace in a 12-dimensional vector space, the projection can be described as a simple linear map.

Using a reduced version of $underline{beta}=begin{pmatrix}beta_{11}beta_{22}beta_{33}end{pmatrix}$ (in the projection step $underline{underline{B}}$ will be flattened to a vector to avoid working with 3D tensors), the projection $underline{beta}'=operatorname{proj}_C(underline{beta})$ can be represented as

$underline{beta}' = underline{underline{P}}cdotunderline{beta}~~~~text{with}~~~~ underline{underline{P}}=frac{1}{3}begin{pmatrix}2&-1&-1-1&2&-1-1&-1&2end{pmatrix}$

The vector $underline{beta}$ and projection map $underline{underline{P}}$ can then be arbitrarily augmented with the remaining entries of $underline{underline{B}}$, as they will not be affected by the projection.

With the above method, my constrained gradient descent is only slightly more expensive than a usual gradient descent method by computing one extra linear operation for each iteration.

The results I have obtained so far with the method seem promising