振幅と位相を合わせて、アンテナゲインを複素数として解きます。ポイントは2点：

「振幅」と「位相」でなく複素数、すなわち「実部と虚部」として解きます。ビジビリティの実部・虚部はそれぞれ正規分布すると期待されるので、最小二乗法にバイアスが乗るのを防ぎます（振幅はRice分布するので最小二乗法にバイアスが乗る）。
refantの位相は0として一般性を失いません。つまりrefantのゲイン$\boldsymbol{G}_0$の虚部は0に固定します。

「振幅を解く」の時と同様に、振幅は適当なスケーリングを援用して

\hat{V}_k = \boldsymbol{G}_j \boldsymbol{G}^*_i \tag{6.1}

を$\boldsymbol{G}$について解くことにします。未知数のベクトルを、$\boldsymbol{G} = (reG_0, reG_1, reG_2, \dots, reG_{Na-1}, imG_1, imG_2, \dots, imG_{Na-1})$ というように、$2N_a - 1$個の実数で表します。$\boldsymbol{G}=R + iIというように$の実部と虚部をそれぞれ$R, I$と書くことにしましょう。観測方程式は行列$P$に未知数を含む非線形方程式ですから、反復法で解きます。初期値を$\boldsymbol{G^0}$とすると残差ベクトルは$r_k = \hat{V}_k - G^0_j G_i^{*0}$で、それに対する修正量ベクトル$\boldsymbol{c}$は$\boldsymbol{r} = P\boldsymbol{c}$を満たし、$\boldsymbol{c} = (P^T P)^{-1} P^T\boldsymbol{r}$によって得られます。行列$P$の成分を書き出すと、

P = \left( \begin{array}{ccccccccc}
 R_1 &  R_0 & 0 & 0 & \cdots &  I_0 & 0 & 0 & \cdots \\
 R_2 &  0 & R_0 & 0 & \cdots &  0 & I_0 & 0 & \cdots \\
 0 & R_2 & R_1 & 0 & \cdots &  I_2 & I_0 & 0 & \cdots \\
 R_3 & 0 & 0 & R_0 & \cdots &  0 & 0 & I_0 & \cdots \\
 0 & R_3 & 0 & R_1 & \cdots &  I_3 & 0 & I_1 & \cdots \\
 0 & 0 & R_3 & R_2 & \cdots &  0 & I_3 & I_2 & \cdots \\
\vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \ddots \\
 I_1 &  -I_0 & 0 & 0 & \cdots &  R_0 & 0 & 0 & \cdots \\
 I_2 &  0 & -I_0 & 0 & \cdots &  0 & R_0 & 0 & \cdots \\
 0 & I_2 & -I_1 & 0 & \cdots &  -R_2 & R_0 & 0 & \cdots \\
 I_3 & 0 & 0 & -I_0 & \cdots &  0 & 0 & R_0 & \cdots \\
 0 & I_3 & 0 & -I_1 & \cdots &  -R_3 & 0 & R_1 & \cdots \\
 0 & 0 & I_3 & -I_2 & \cdots &  0 & -R_3 & R_2 & \cdots \\
\vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \ddots \\
\end{array} \right) = \left( \begin{array}{cc}
A & B \\
C & D \\
\end{array} \right) \tag{6.2}

となり、4つの部分行列$A, B, C, D$で表されます。

A =  \left( \begin{array}{ccccc}
 R_1 &  R_0 & 0 & 0 & \cdots \\
 R_2 &  0 & R_0 & 0 & \cdots \\
 0 & R_2 & R_1 & 0 & \cdots  \\
 R_3 & 0 & 0 & R_0 & \cdots \\
 0 & R_3 & 0 & R_1 & \cdots  \\
 0 & 0 & R_3 & R_2 & \cdots  \\
\vdots & \vdots & \vdots & \vdots & \ddots  \\
\end{array} \right) \\
B = \left( \begin{array}{cccc}
  I_0 & 0 & 0 & \cdots \\
  0 & I_0 & 0 & \cdots \\
 I_2 & I_0 & 0 & \cdots \\
  0 & 0 & I_0 & \cdots \\
 I_3 & 0 & I_1 & \cdots \\
 0 & I_3 & I_2 & \cdots \\
\vdots & \vdots & \vdots & \ddots \\
\end{array} \right) \\
C = \left( \begin{array}{ccccc}
 I_1 &  -I_0 & 0 & 0 & \cdots  \\
 I_2 &  0 & -I_0 & 0 & \cdots  \\
 0 & I_2 & -I_1 & 0 & \cdots  \\
 I_3 & 0 & 0 & -I_0 & \cdots \\
 0 & I_3 & 0 & -I_1 & \cdots \\
 0 & 0 & I_3 & -I_2 & \cdots \\
\vdots & \vdots & \vdots & \vdots & \ddots \\
\end{array} \right) \\
D = \left( \begin{array}{cccc}
R_0 & 0 & 0 & \cdots \\
0 & R_0 & 0 & \cdots \\
-R_2 & R_0 & 0 & \cdots \\
0 & 0 & R_0 & \cdots \\
-R_3 & 0 & R_1 & \cdots \\
0 & -R_3 & R_2 & \cdots \\
\vdots & \vdots & \vdots & \ddots \\
\end{array} \right) \\

行列計算の省力化

部分行列で表すと以下のようになります。

P^T P = \left( \begin{array}{cc}
A^T A + C^TC & A^T B + C^T D \\
B^T A + D^T C & B^TB + D^T D
\end{array} \right) \tag{6.3}

$A^TA$の計算は、「振幅を解く」の式5.7でやっていますのでそのまま使います。$C^TC$は$A^TA$とよく似ていて、

C^TC = \left( \begin{array}{ccccc}
\sum_k I^2_k - I^2_0 & -I_0 I_1 & -I_0 I_2 & -I_0 I_3 & \cdots \\
-I_0 I_1 & \sum_k I^2_k - I^2_1 &  -I_1 I_2 & -I_1 I_3 & \cdots \\
-I_0 I_2 & -I_1 I_2 & \sum_k I^2_k - I^2_2 & -I_2 I_3 & \cdots \\
-I_0 I_3 & -I_1 I_3 & -I_2 I_3 & \sum_k I^2_k - I^2_3 & \cdots \\
\vdots & \vdots & \vdots & \vdots &\ddots
\end{array} \right) \tag{6.4}

となります。$B^TA$, $D^TC$の計算はそれぞれ以下の通りです。

B^TA = \left( \begin{array}{ccccc}
R_1 I_0 & \boldsymbol{R}\cdot\boldsymbol{I} - R_1 I_1 & R_1 I_2 & R_1 I_3 & \cdots \\
R_2 I_0 & R_2 I_1  & \boldsymbol{R}\cdot\boldsymbol{I} - R_2 I_2 &  R_2 I_3 & \cdots \\
R_3 I_0 & R_3 I_1 & R_3 I_2 & \boldsymbol{R}\cdot\boldsymbol{I} - R_3 I_3 & \cdots \\
\vdots & \vdots & \vdots & \vdots &\ddots
\end{array} \right) \tag{6.5}

D^TC = \left( \begin{array}{ccccc}
R_0 I_1 & -\boldsymbol{R}\cdot\boldsymbol{I} + R_1 I_1 & R_2 I_1 & R_3 I_1 & \cdots \\
R_0 I_2 & R_1 I_2  & -\boldsymbol{R}\cdot\boldsymbol{I} + R_2 I_2 &  R_3 I_2 & \cdots \\
R_0 I_3 & R_1 I_3 & R_2 I_3 & -\boldsymbol{R}\cdot\boldsymbol{I} + R_3 I_3 & \cdots \\
\vdots & \vdots & \vdots & \vdots &\ddots
\end{array} \right) \tag{6.6}

なお、$I_0 = 0$です。また、$A^TB = (B^TA)^T$, $C^TD = (D^TC)^T$によって求まります。以上をまとめて、$P^TP$を得る関数をPythonで書くと、以下のようになります。

def ATAmatrix(Gain):  # Gain is a vector of antenna-based gain amplitude (real)
    antNum = len(Gain); normG = Gain.dot(Gain)
    PTP = np.zeros([antNum, antNum]) + Gain
    for ant_index in range(antNum):
        PTP[ant_index,:] *= Gain[ant_index]
        PTP[ant_index, ant_index] = normG - Gain[ant_index]**2
    #
    return PTP
#
def CTCmatrix(Gain):  # Gain is a vector of antenna-based gain amplitude (imag)
    antNum = len(Gain); normG = Gain.dot(Gain)
    PTP = np.zeros([antNum, antNum]) + Gain
    for ant_index in range(antNum):
        PTP[ant_index,:] *= (-Gain[ant_index])
        PTP[ant_index, ant_index] = normG - Gain[ant_index]**2
    #
    return PTP
#
def ATBmatrix(Gain):  # Gain is a vector of antenna-based complex gain
    antNum = len(Gain); normG = Gain.real.dot(Gain.imag)
    PTP = np.zeros([antNum, antNum]) + Gain.imag
    for ant_index in range(antNum):
        PTP[ant_index,:] = Gain.real[ant_index]* Gain.imag + Gain.imag[ant_index]* Gain.real
        PTP[ant_index, ant_index] = 0.0
    #
    return PTP[1:antNum]
#
def PMatrix(CompSol):
    antNum = len(CompSol); matSize = 2*antNum-1
    PM = np.zeros([matSize, matSize])
    PM[0:antNum][:,0:antNum]   = ATAmatrix(CompSol.real) + CTCmatrix(CompSol.imag) # ATA + CTC
    PM[antNum:matSize][:,antNum:matSize] = ATAmatrix(CompSol.imag)[1:antNum][:,1:antNum] + CTCmatrix(CompSol.real)[1:antNum][:,1:antNum] # BTB + DTD
    PM[antNum:matSize][:,0:antNum]   = ATBmatrix(CompSol) # ATB + CTD
    PM[0:antNum][:,antNum:matSize]   = PM[antNum:matSize][:,0:antNum].T # BTA + DTC
    return PM
#

引数のCompSolはアンテナゲイン（複素数）初期値のベクトルです。返り値の行列PMは成分が実数でサイズが$(2N_a - 1, 2N_a - 1)$です。

次に$P^T r$ を計算します。以下のPython関数PTdotR()は、引数にアンテナゲイン（複素数）初期値のベクトルCompSolと残差ベクトル (複素数) Cresid をとり、$P^T r$を計算して返します。

def PTdotR(CompSol, Cresid):
    antNum = len(CompSol)
    blNum = antNum* (antNum-1) / 2
    ant0, ant1= np.array(ANT0[0:blNum]), np.array(ANT1[0:blNum])
    PTR = np.zeros(2*antNum)
    for ant_index in range(antNum):
        index0 = np.where(ant0 == ant_index)[0].tolist()
        index1 = np.where(ant1 == ant_index)[0].tolist()
        PTR[range(ant_index)]          += (CompSol[ant_index].real* Cresid[index0].real + CompSol[ant_index].imag* Cresid[index0].imag)
        PTR[range(ant_index+1,antNum)] += (CompSol[ant_index].real* Cresid[index1].real - CompSol[ant_index].imag* Cresid[index1].imag)
        PTR[range(antNum, antNum+ant_index)]    += (CompSol[ant_index].imag* Cresid[index0].real - CompSol[ant_index].real* Cresid[index0].imag)
        PTR[range(antNum+ant_index+1,2*antNum)] += (CompSol[ant_index].imag* Cresid[index1].real + CompSol[ant_index].real* Cresid[index1].imag)
    #
    return PTR[range(antNum) + range(antNum+1, 2*antNum)]
#

まとめのコード

上記のPMatrix()およびPTdotR()を使って、反復法でアンテナ複素ゲインを解く関数 gainComplex() を下記に示します。引数のbl_vis は基線ベースのビジビリティ（複素数）で、Canonical orderingに従っているものとします。niterは反復回数ですが、2回程度で十分に収束するでしょう。$P^{T}P$は実対称行列ですので、Cholesky分解によって下三角行列Lを得て方程式を解くことにします。

def gainComplex( bl_vis, niter=2 ):
    blNum  =  len(bl_vis)
    antNum =  Bl2Ant(blNum)[0]
    ant0, ant1, kernelBL = ANT0[0:blNum], ANT1[0:blNum], KERNEL_BL[range(antNum-1)].tolist()
    CompSol = np.zeros(antNum, dtype=complex)
    #---- Initial solution
    CompSol[0] = sqrt(abs(bl_vis[0])) + 0j
    CompSol[1:antNum] = bl_vis[kernelBL] / CompSol[0]
    #----  Iteration
    for iter_index in range(niter):
        PTP        = PMatrix(CompSol)
        L          = np.linalg.cholesky(PTP)         # Cholesky decomposition
        Cresid     = bl_vis - CompSol[ant0]* CompSol[ant1].conjugate()
        t          = np.linalg.solve(L, PTdotR(CompSol, Cresid))
        correction = np.linalg.solve(L.T, t)
        CompSol    = CompSol + correction[range(antNum)] + 1.0j* np.append(0, correction[range(antNum, 2*antNum-1)])
    #
    return CompSol
#

以上です。

もくじにもどる