干渉計観測の遅延を解く

#はじめに
遅延$\tau$は位相と角周波数$2\pi \nu$との間に$\phi = \phi_0 + 2 \pi \nu \tau$という関係を持ちます。基線毎のビジビリティで得られる遅延 $ \hat{\tau}$ は、アンテナ毎の遅延 $ \tau $ に対して、

```math
\hat{\tau}_k = \tau_j - \tau_i \tag{3.1}
```

で表されます。$i, j, k$は[canonical ordering](http://qiita.com/kamenoseiji/items/ed6fd86641d088da9fff) に従うものとします。この関係は全てのアンテナに共通の遅延オフセットを与えても成り立ちますから、基準となるアンテナ (refant) の遅延を0と置いても一般性を失いません。

##最小二乗法
基線毎に測定された遅延$\hat{\tau}_k$からアンテナ毎の遅延$\tau_i$を最小二乗法で解く方程式は

```math
\left(
\begin{array}{c}
\hat{\tau}_0  \\
\hat{\tau}_1  \\
\hat{\tau}_2  \\
\hat{\tau}_3  \\
\hat{\tau}_4  \\
\hat{\tau}_5  \\
\vdots
\end{array}
\right) = \left( \begin{array}{cccc}
 1 & 0 & 0 & \cdots \\
 0 & 1 & 0 & \cdots \\
-1 & 1 & 0 & \cdots \\
 0 & 0 & 1 & \cdots \\
-1 & 0 & 1 & \cdots \\
 0 & -1 & 1 & \cdots \\
\vdots & \vdots & \vdots & \ddots
\end{array} \right)
\left(
\begin{array}{c}
\tau_1  \\
\tau_2  \\
\tau_3  \\
\vdots
\end{array}
\right) \tag{3.2}
```

という具合になります。$\tau_0=0$と置いたので未知数は$N_a - 1$個となります。真ん中の$(N_a - 1) \times N_b$サイズの行列を$P$とおき、観測量$\hat{\tau}_k$のベクトルを$\boldsymbol{Y}$、未知数$\tau_i$のベクトルを$\boldsymbol{X}$とおけば、アンテナ毎の遅延を最小二乗の意味で得る連立方程式は

```math
(P^TP) \boldsymbol{X} = P^T \boldsymbol{Y} \tag{3.3}
```

となります。従って$\boldsymbol{X} = (P^TP)^{-1} P^T \boldsymbol{Y}$ですから$\boldsymbol{X}$は求められます。

## 計算の省力化
アンテナ数・基線数が大きくなると$P$の演算量は膨大になります。ALMAではアンテナ数が最大64, 基線数が最大2016ですので、これを解くのは大変です。幸い、$P^TP$の計算はcanonical orderingを採用していると以下のようにシンプルになります。

```math
P^TP = \left( \begin{array}{ccccc}
 N_a - 1 & -1 & -1 & -1 & \cdots \\
-1 & N_a - 1 & -1 & -1 & \cdots \\
-1 & -1 & N_a-1 & -1 & \cdots \\
-1 & -1 & -1 & N_a-1 & \cdots \\
\vdots & \vdots & \vdots & \vdots &\ddots
\end{array} \right) \tag{3.4}
```

要するに対角成分が$N_a-1$で、それ以外の成分は$-1$です。これは実対称行列ですからCholesky分解を用いて$LL^T = P$と表せる下三角行列$L$を得て解けます。さらに言うと$(P^TP)$の逆行列は簡単に直接計算できて、

```math
(P^TP)^{-1}= \frac{1}{N_a}\left( \begin{array}{ccccc}
 2 & 1 & 1 & 1 & \cdots \\
 1 & 2 & 1 & 1 & \cdots \\
 1 & 1 & 2 & 1 & \cdots \\
 1 & 1 & 1 & 2 & \cdots \\
\vdots & \vdots & \vdots & \vdots &\ddots
\end{array} \right) \tag{3.5}
```

となります。要するに対角成分が$\frac{2}{N_a}$, それ以外の成分が$\frac{1}{N_a}$となる$(N_a - 1) \times (N_a - 1)$サイズの行列です。Pythonコードで書くとantNumをアンテナ数として実質1行で済みます。

```Python
import numpy as np
PTP_inv = (np.diag(np.ones(antNum - 1)) + 1.0) / antNum
```

残るは$P^T \boldsymbol{Y}$の計算です。巨大な疎行列になる$P^T$を計算するのは避けるにしても、基線数のループは避けたいので、以下のようにアンテナ数のループ内でベクトル演算をします。

```Python
ant0, ant1 = np.array(ANT0[0:blNum]), np.array(ANT1[0:blNum])
PTY = np.zeros(antNum - 1)
for ant_index in range(1, antNum):
  index0 = np.where(ant0 == ant_index)[0].tolist()
  index1 = np.where(ant1 == ant_index)[0].tolist()
  PTY[ant_index - 1] += np.sum(Y[index0]) - np.sum(Y[index1])
#
```

なお、リストANT0, ANT1はあらかじめこの記事[Canonical Orderingによる干渉計の基線番号付け](http://qiita.com/kamenoseiji/items/ed6fd86641d088da9fff)を参照して作っておきます。index0は「ant_indexを基線終点にもつ
このようにしてできた$(P^TP)^{-1}$と$P^TY$の内積をとれば、未知数のベクトル$X$が得られます。

```Python
Delay_ant = [0.0] + (PTP_inv.dot(PTY)).tolist()
```

リストの頭に0.0を挿入しているのはrefantを含めたリストにしているからです。このようにしておくと、基線ごとの遅延を

```Python
Delay_bl = Delay_ant[ANT1] - Delay_ant[ANT0]
```
のようにして得られます。

##まとめのコード
基線ベースの遅延からアンテナベースの遅延を得るPythonの関数を以下に記します。引数のbl_delayは基線ベースの遅延を並べたnumpy arrayで、[Canonical Ordering](http://qiita.com/kamenoseiji/items/ed6fd86641d088da9fff) で並んでいる必要があります。返り値はアンテナベースの遅延で、先頭はrefantですので0です。

```Python
def cldelay_solve(bl_delay):
    blNum = len(bl_delay); antNum = Bl2Ant(blNum)[0]
    ant0, ant1 = np.array(ANT0[0:blNum]), np.array(ANT1[0:blNum])
    PTP_inv = (np.diag(np.ones(antNum - 1)) + 1.0) / antNum
    PTY = np.zeros(antNum - 1)
    for ant_index in range(1, antNum):
        index0 = np.where(ant0 == ant_index)[0].tolist()
        index1 = np.where(ant1 == ant_index)[0].tolist()
        PTY[ant_index - 1] += np.sum(bl_delay[index0]) - np.sum(bl_delay[index1])
    #
    return np.array( [0.0] + (PTP_inv.dot(PTY)).tolist() )
#
```

以上です。

[もくじにもどる](http://qiita.com/kamenoseiji/items/0862bf3b3b57d203009a)