振幅と位相を合わせて、アンテナゲインを複素数として解きます。ポイントは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
#
以上です。