[PYTHON] Solving the complex gain of interferometer observations

Match the amplitude and phase and solve the antenna gain as a complex number. 2 points:

1. Complex numbers, not "amplitude" and "phase" Solve as "part and imaginary part". Since the real and imaginary parts of the visibility are expected to be normally distributed, it prevents the least squares method from being biased (amplitude is a Rice distribution. Therefore, the least squares method is biased).
2. The phase of refant is set to 0 and does not lose generality. In other words, the imaginary part of the gain $ \ boldsymbol {G} _0 $ of the refund is fixed at 0.

As in the case of "Solving the amplitude", the amplitude uses appropriate scaling.

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

Let's solve for $ \ boldsymbol {G} $. The unknown vector is $ 2N_a, such as $ \ boldsymbol {G} = (reG_0, reG_1, reG_2, \ dots, reG_ {Na-1}, imG_1, imG_2, \ dots, imG_ {Na-1}) $. --Represented by 1 $ real number. Let's write the real and imaginary parts of $ as $ R and I $, respectively, like $ \ boldsymbol {G} = R + iI. Since the observation equation is a nonlinear equation containing unknowns in the matrix $ P $, it is solved by the iterative method. If the initial value is $ \ boldsymbol {G ^ 0} $, the residual vector is $ r_k = \ hat {V} _k --G ^ 0_j G_i ^ {* 0} $, and the correction amount vector for it $ \ boldsymbol {c } $ Satisfies $ \ boldsymbol {r} = P \ boldsymbol {c} $ and is obtained by $ \ boldsymbol {c} = (P ^ TP) ^ {-1} P ^ T \ boldsymbol {r} $ .. When you write out the components of the matrix $ 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}

And is represented by four submatrix $ 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) \\

Labor saving in matrix calculation

It is expressed as a submatrix as follows.

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}

The calculation of $ A ^ TA $ is done by the formula 5.7 of "Solving the amplitude", so use it as it is. $ C ^ TC $ is very similar to $ 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}

It will be. The calculations for $ B ^ TA $ and $ D ^ TC $ are as follows.

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}

Note that $ I_0 = 0 $. It is also calculated by $ A ^ TB = (B ^ TA) ^ T $, $ C ^ TD = (D ^ TC) ^ T $. Summarizing the above, if you write a function to get $ P ^ TP $ in Python, it will be as follows.

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
#

The argument CompSol is a vector of the initial antenna gain (complex number). The return matrix PM has real components and a size of $ (2N_a ―― 1, 2N_a ―― 1) $.

Then calculate $ P ^ T r $. The following Python function PTdotR () takes the antenna gain (complex number) initial value vector CompSol and the residual vector (complex number) Cresid as arguments, and calculates and returns $ 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)]
#

Summary code

The function gainComplex () that solves the antenna complex gain by iterative method using PMatrix () and PTdotR () above is shown below. The argument bl_vis is a baseline-based visibility (complex number) and is assumed to follow Canonical ordering. The niter is the number of iterations, but about 2 will be enough to converge. Since $ P ^ {T} P $ is a real symmetric matrix, we will solve the equation by obtaining the lower triangular matrix L by Cholesky decomposition.

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
#

that's all.

Return to the lottery