# Jaccobi Eigenvalue Algorithm

應用時機: 求解**實對稱矩陣**的特徵值與特徵向量，在矩陣階數小於10以下，相對於QR分解快速。

## 原理
如果矩陣 $A$ 是實對稱矩陣，可以透過一系列精心設計過的旋轉矩陣組成正交矩陣 $V$，使得 $D=V^TAV$ 為對角矩陣。該變換實則為相似變換，所以矩陣 $D$ 的對角元素為 $A$ 的特徵值。

Jacobi特徵法由 $D = R_k^T...R_1^TAR_1...R_k$ 形式組成。每個正交變換(**Jacobi rotations**)都是在旋轉平面，主要是在消滅**非對角**元素。

最後該方法收斂的證明詳見 [線代啟示錄-Jacobi 特徵值算法](https://ccjou.wordpress.com/2015/09/17/jacobi-%E7%89%B9%E5%BE%B5%E5%80%BC%E7%AE%97%E6%B3%95/)。

:::spoiler Given Rotation 如何消滅非對角元素?
詳見[Givens 旋轉於 QR 分解的應用](https://ccjou.wordpress.com/2010/02/18/givens-%E6%97%8B%E8%BD%89%E6%96%BC-qr-%E5%88%86%E8%A7%A3%E7%9A%84%E6%87%89%E7%94%A8/)

舉2維空間作為例子，令$c=cos\theta$,$s=sin\theta$。Given Rotation旨在以$x_p$為旋轉軸，消滅$x_q$，如下所述：

$$
\begin{bmatrix}
c & -s\\
s & c \\
\end{bmatrix}
\begin{bmatrix}
x_p \\
x_q \\
\end{bmatrix}
\begin{bmatrix}
r \\
0 \\
\end{bmatrix}
$$

:::

:::spoiler Jacobi Rotation 如何消滅非對角元素?
根據上述式子 $D=V^TAV$，分成兩個部分介紹：$D$ 對角矩陣 以及 $V$ 正交矩陣。
==**(1) $D$ 對角矩陣由來**==

舉2維空間作為例子，精心設計$\theta$，使得$d_{pq}=d_{qp}=0$，求$\theta$：

$$
D=
\begin{bmatrix}
d_{pp} & d_{pq} \\
d_{qp} & d_{qq} \\
\end{bmatrix}
=R^TAR=
\begin{bmatrix}
cos\theta & sin\theta \\
-sin\theta & cos\theta \\
\end{bmatrix}
\begin{bmatrix}
a_{pp} & a_{pq} \\
a_{pq} & a_{qq} \\
\end{bmatrix}
\begin{bmatrix}
cos\theta & -sin\theta \\
sin\theta & cos\theta \\
\end{bmatrix}
$$

乘開可得，並令其為0：

$$
d_{pq}=cos\theta sin\theta (a_{qq}-a_{pp}) + (cos^2\theta - sin^2\theta) a_{pq}=d_{qp}=0
$$

透過倍角公式 $cos2\theta = cos^2\theta - sin^2\theta$，$sin2\theta = 2cos\theta sin\theta$ 以及 $tan2\theta = \frac{2tan\theta}{1-tan^2\theta}$ ，可得：

$$
tan2\theta = \frac{2a_{pq}}{a_{pp}-a_{qq}}=\frac{2tan\theta}{1-tan^2\theta}
$$

令其 ==$\phi=1/tan2\theta$== 以及 ==$t=tan\theta$==，

$$
\phi=\frac{ 1 - t^2 }{ 2t }=\frac{ a_{pp} - a_{qq} }{ 2a_{pq} } \implies t^2+2t\phi -1 = 0
$$

選擇絕對值較小的解 $t_{min}=sgn(\phi)(-|\phi|+\sqrt{\phi^2+1})$，為了避免兩個相近大小的數相減造成[catastrophic cancellation](https://zh.wikipedia.org/wiki/%E7%81%BE%E9%9A%BE%E6%80%A7%E6%8A%B5%E6%B6%88)，對 $t_{min}$ 乘上 $\frac{|\phi+\sqrt{\phi^2+1}|}{|\phi+\sqrt{\phi^2+1}|}$ ，最後得出：

<center>

==$t=\frac{sgn(\phi)}{|\phi|+\sqrt{\phi^2+1}}$==

</center>

但如果 $\phi >> 1$ ，恐會讓 $\phi^2$ 造成overflow，因此在該狀況下讓 ==$t=\frac{1}{2\phi}$==。

另外Jacobi特徵法的參數值如下，
<center>

==$c=cos\theta=\frac{1}{\sqrt{1+t^2}}$== , ==$s=sin\theta=t \cdot c$==

</center>

當$2\theta = arctan\frac{2a_{pq}}{a_{qq}-a_{pp}}$時，使得$d_{pq}=d_{qp}=0$。證明了透過jacobi rotation二階實對稱矩陣可正交對角化。

**擴展至n維狀況**:
正交相似變換定義為$A'=R_{pq}^T \cdot A \cdot R_{pq}$，其中旋轉矩陣如下所示：

$$
R_{pq}=
\begin{bmatrix}
1 & & & & & & \\
&...& & & & & \\
& & c &...&-s& & \\
& & : & 1 & :& & \\
& & s &...& c& & \\
& & & & &... & \\
& & & & & & 1\\
\end{bmatrix}
$$

上述操作可知 $R_{pq}^T \cdot A$ **改變**了 $A$ 的==p,q列== 以及$A \cdot R_{pq}$則**改變**了 $A$的==p,q行==。 如下圖所示，黃線處為變動項：

<center>

<img width=45% src=https://hackmd.io/_uploads/HkrP3wEQj.png>

</center>

乘開後便有以下規律:

<center>

\begin{cases}
a'_{rp}=ca_{rp}-sa_{rq}, & \text{$r \neq p, r\neq q$} \\
a'_{rq}=ca_{rq}+sa_{rp}, & \text{$r \neq p, r\neq q$} \\
a'_{pp}=c^2a_{pp}+s^2a_{qq}-2sc\cdot a_{pq} \\
a'_{qq}=s^2a_{pp}+c^2a_{qq}+2sc\cdot a_{pq} \\
a'_{pq} = (c^2-s^2)a_{pq} + sc(a_{pp}-a_{qq})=0
\end{cases}

</center>

整理成:

$$
\begin{bmatrix}
a'_{rp}\\
a'_{rq}\\
\end{bmatrix}=
\begin{bmatrix}
c & -s\\
s & c\\
\end{bmatrix}
\begin{bmatrix}
a_{rp}\\
a_{rq}\\
\end{bmatrix},\quad r \neq p, r\neq q \\
a'_{pp}=a_{pp}-t \cdot a_{pq} \\
a'_{qq}=a_{qq}+t \cdot a_{pq} \\
a'_{pq}=0
$$

==**(2) $V$ 正交矩陣由來**==
從上述迭代式，可知正交矩陣 $V$ 之形式：
$$
A=VDV^T \Longleftrightarrow D=V^TAV=(R_k^T \cdots R_1^T)A(R_1 \cdots R_k)
$$
代表 $V=I\cdot R_1\cdots R_k$，公式整理如下：

$$
\begin{bmatrix}
v'_{rp}\\
v'_{rq}\\
\end{bmatrix}=
\begin{bmatrix}
c & -s\\
s & c\\
\end{bmatrix}
\begin{bmatrix}
v_{rp}\\
v_{rq}\\
\end{bmatrix}
$$
:::

::: spoiler 那些Jacobi Rotations以C語言實現的坑

因為$A$矩陣為實對稱矩陣，實質上也只需要尋訪上三角$n(n-1)/2$個。而因為ｃ語言實現，故旋轉就分成三種case。

<center>
<img width=45% src=https://hackmd.io/_uploads/HJxTlgBXi.png>
</center>

//1. Case of rotations 0 ≤ j< p.
for (int j = 0; j<ip; j++)
    Rotate(ref A, j, ip, j, iq, c, s);

//2. Case of rotations p+1 ≤ j< q.
for (int j = ip+1; j<iq; j++)
    Rotate(ref A, ip, j, j, iq, c, s);

//3. Case of rotations q+1 ≤ j< n.
for (int j = iq+1; j<n; j++)
    Rotate(ref A, ip, j, iq, j, c, s);

void Rotate(double[,] A, int r1, int c1, int r2, int c2, double c, double s)
{
    double num = A[r1, c1];
    double num2 = A[r2, c2];
    A[r1, c1] = num * c - num2 * s;
    A[r2, c2] = num * s + num2 * c;
}

:::

::: spoiler Cyclic Jacobi method 代碼

//From Numerical Recipes in C. Second Edition
public static void JacobiEigenvalueAlgorithm(
    double[,] A,
    ref double[] W,
    ref double[,] V,
    out int totalitera)
{
    int n = A.GetLength(0);

    totalitera = 0;
    int Maxitera = n * n * 30;

    for (int i = 0; i < n; i++)
    {
        //Initialize V to the identity matrix
        if (V != null)
            V[i, i] = 1.0;
        //Initialize W to the diagonal of A matrix
        W[i] = A[i, i];
    }

    double Magnitude_OffDiag, h, t, c, s, diagAddtemr;
    double threshold;
    double sum_OffDiag;
    for (int i = 0; i<Maxitera; i++)
    {
        //Sum off-diagonal elements.
        sum_OffDiag=0;
        for (int ip = 0; ip<n-1; ip++)
            for (int iq = ip+1; iq<n; iq++)
                sum_OffDiag +=Math.Abs(A[ip, iq]);

        // if the sum of off-diagonal elements is zeros,
        //  then it reaches convergence.
        if (sum_OffDiag <2.2204460492503131E-16)
            break;

        threshold = (i<4) ? 0.2*sum_OffDiag/(n*n) : 0;
        for (int ip = 0; ip<n-1; ip++)
        {
            for (int iq = ip+1; iq<n; iq++)
            {
                Magnitude_OffDiag=100.0*Math.Abs(A[ip, iq]);
                //After four sweeps,
                //skip the rotation 
                //if the off-diagonal element is small 
                //(W[ip] and W[iq] >>  100*|A[ip,iq]| ).
                if (i>3
                    && (Math.Abs(W[ip])+Magnitude_OffDiag)==Math.Abs(W[ip])
                    && (Math.Abs(W[iq])+Magnitude_OffDiag)==Math.Abs(W[iq])
                    )
                {
                    A[ip, iq]=0;
                }
                else if (Math.Abs(A[ip, iq])>threshold)
                {
                    h=W[iq]-W[ip];
                    if ((Math.Abs(h)+Magnitude_OffDiag)==Math.Abs(h))
                        t=A[ip, iq]/h;
                    else
                    {
                        //1. phi = (A[p,p] -A[q,q]) / (2*A[p,q]);
                        //2. t = tan(theta) =  sgn(phi) /(|phi|+sqrt(phi^2 +1));
                        double phi = 0.5 * h/A[ip, iq];
                        t = 1.0/(Math.Abs(phi)+Math.Sqrt(1.0+phi*phi));
                        if (phi<0)
                            t=-t;
                    }
                    //3. c= cos(theta) = 1 / sqrt(1 + tan(theta)^2 )
                    //4. s= sin(theta) = tan(theta) * cos(theta)
                    c=1.0/Math.Sqrt(1+t*t);
                    s=t*c;

                    // Follow eq(11.1.14)
                    diagAddtemr=t*A[ip, iq];
                    W[ip] -= diagAddtemr;
                    W[iq] += diagAddtemr;

                    A[ip, iq]=0;
                    //Case of rotations 0 ≤ j< p.
                    for (int j = 0; j<ip; j++)
                        Rotate(ref A, j, ip, j, iq, c, s);

                    //Case of rotations p+1 ≤ j< q.
                    for (int j = ip+1; j<iq; j++)
                        Rotate(ref A, ip, j, j, iq, c, s);

                    //Case of rotations q+1 ≤ j< n.
                    for (int j = iq+1; j<n; j++)
                        Rotate(ref A, ip, j, iq, j, c, s);

                    if (V != null)
                    {
                        for (int j = 0; j<n; j++)
                            Rotate(ref V, j, ip, j, iq, c, s);
                    }

                    totalitera++;
                }
            }
        }
    }

    //sorting by descent via bubble sort
    for (int j = 0; j < n - 1; j++)
    {
        int key = j;
        for (int i = j + 1; i < n; i++)
        {
            if (Math.Abs(W[key]) < Math.Abs(W[i]))
                key = i;

        }

        if (j != key)
        {
            swap(ref W[key], ref W[j]);
            if (V!=null)
            {
                for (int i = 0; i < n; i++)
                    swap(ref V[i, key], ref V[i, j]);
            }
        }
    }
}
private static void Rotate(
    ref double[,] A,
    int row1, int col1,
    int row2, int col2,
    double c, double s)
{
    double num = A[row1, col1];
    double num2 = A[row2, col2];
    A[row1, col1] = num * c - num2 * s;
    A[row2, col2] = num * s + num2 * c;
}

:::

## Reference:

1. Numerical Recipes in C , Second Edition (1992), Ch11.
2. Jacobi eigenvalue algorithm, WIKI
3. [線代啟示錄-Jacobi 特徵值算法](https://ccjou.wordpress.com/2015/09/17/jacobi-%E7%89%B9%E5%BE%B5%E5%80%BC%E7%AE%97%E6%B3%95/)
4. [Givens 旋轉於 QR 分解的應用](https://ccjou.wordpress.com/2010/02/18/givens-%E6%97%8B%E8%BD%89%E6%96%BC-qr-%E5%88%86%E8%A7%A3%E7%9A%84%E6%87%89%E7%94%A8/)

###### tags: `Math`