bilinear complex least squares method I want to minimise S, which is given by,
$S = \sum_{i,j} \frac{1}{\sigma^2_{ij}}|V_{ij} - g_i g^*_j V^{model}_{ij}|^2$,
where $V_{ij}$ and $V^{model}_{ij}$ are the observed and model complex visibilities (this is a (M x 1) vector) both of which are known, $\sigma_{ij}$ is the error and $g_i$ and $g_j$ are the complex gains that I am minimising S for.
In my case I am only interested to solve for the phases of the complex gains assuming their amplitude is 1. I can therefore write the following:
$V^{model, corrupted}_{ij} = g_i g^*_j V^{model}_{ij} = e^{i (\phi_i - \phi_j)} V^{model}_{ij} = e^{i f \phi} V^{model}_{ij}$,
where f is a (M x N) matrix and $\phi$ is a (N x 1) vector of the gain phases (for N=4 antennas I have M=N(N-1)/2 visibilities), so that
$\begin{pmatrix} 
\phi_1 - \phi_2 \\ 
\phi_1 - \phi_3 \\ 
\phi_1 - \phi_4 \\ 
\phi_2 - \phi_3 \\
\phi_2 - \phi_4 \\
\phi_3 - \phi_4 \\
\end{pmatrix} = 
\begin{pmatrix}
    1 & -1 & 0 & 0 \\
    1 & 0 & -1 & 0 \\
    1 & 0 & 0 & -1 \\
    0 & 1 & -1 & 0 \\
    0 & 1 & 0 & -1 \\
    0 & 0 & 1 & -1 \\
\end{pmatrix} \cdot
\begin{pmatrix} \phi_1 & \phi_2 & \phi_3 & \phi_4 \end{pmatrix} = f\phi$
I can now minimise the expression S, using for example a gradient descent algorithm, for the vector $\phi$. However, as you can imagine as the number of antennas increases the problem becomes more and more computationally demanding.
I was wondering if there is a more efficient way to approach my problem, either analytic (doubtful) or with a more sophisticated algorithm that takes advantage of the fact that my solution has the form that I want (i.e. $g_j = e^{i \phi_j}$).
 A: I assume that the $\sigma_{ij}$ are known.
Let $\Sigma$ denote the matrix with entries $\sigma_{ij}$.  Let $U$ denote the matrix
$$
U = \exp\operatorname{diag}(\phi_1,\phi_2,\phi_3) = \pmatrix{e^{i\phi_1 }\\ & e^{i \phi_2}\\ & & e^{i \phi_3}}.
$$
Note that the matrix whose entries are $V_{ij}g_ig_j^*$ can be written as $UVU^*$.  With that, your objective function can be expressed as
$$
S = \sum_{i,j} \frac{1}{\sigma^2_{ij}}|V_{ij} - g_i g^*_j V^{model}_{ij}|^2 = \|\Sigma \odot (V - UV^{\text{model}}U^*)\|_F^2,
$$
where $\odot$ denotes a Hadamard product. With that, a gradient with respect to the vector $\phi = (\phi_1,\phi_2,\phi_3)$ can be computed.

We have:
$$
dS = d\|\Sigma\odot(V - UV^MU^*)\|_F^2 = 
2 \operatorname{Re}\operatorname{tr}[(\Sigma\odot(V - UV^MU^*))^*d(\Sigma\odot(V - UV^MU^*))].
\\
d[\Sigma\odot(V - UV^MU^*)] = \Sigma \odot d(UV^MU^*).\\
d [UV^MU^*] = U(\phi + d\phi)V^MU^*(\phi + d\phi) - UV^MU^* \\
= U V^M[\operatorname{diag}(i e^{i\phi}d\phi)]^*
+ [\operatorname{diag}(i e^{i\phi}d\phi)]V^MU^*.\\
$$
To find the entries of the gradient, plug in $dx = e_j$ (the standard basis vectors). We have $dx = e_j \implies$
$$
d [UV^MU^*] = UV^M\operatorname{diag}(i e^{i\phi_j} e_j)^* + 
\operatorname{diag}(i e^{i\phi_j} e_j)V^MU^*\\
= -ie^{-i\phi_j} U (V^M e_j)e_j^T  + 
ie^{i\phi_j}e_j(e_j^T V^M)V^MU^*.\\
d[\Sigma\odot(V - UV^MU^*)] = 
\Sigma \odot [-ie^{-i\phi_j} U (V^M e_j)e_j^T  + 
ie^{i\phi_j}e_j(e_j^T V^M)V^MU^*].\\
dS =\\ 
2 \operatorname{Re}\operatorname{tr}[(\Sigma\odot(V - UV^MU^*))^*(\Sigma \odot [-ie^{-i\phi_j} U (V^M e_j)e_j^T  + 
ie^{i\phi_j}e_je_j^TV^MU^*])]
\\ = 
2 \operatorname{Re}\operatorname{tr}[\Sigma^{\odot2} \odot ((V - UV^MU^*)^*[-ie^{-i\phi_j} U (V^M e_j)e_j^T  + 
ie^{i\phi_j}e_je_j^T V^MU^*])]
\\ = 
2 \operatorname{Re}[
-ie^{-i\phi_j}
e_j^T( \Sigma^{\odot2} \odot [(V - UV^MU^*)^* U V^M ])e_j\\
\quad +ie^{-i\phi_j}
e_j^T( \Sigma^{\odot2} \odot [(V - UV^MU^*)^*  V^MU^*])e_j].
$$
