循環矩陣傅里葉對角化
All circulant matrices are made diagonal by the Discrete Fourier Transform (DFT), regardless of the generating vector x.
任意循環矩陣可以被傅里葉變換矩陣對角化。
文獻中，一般用如下方式表達這一概念：
X
=
C
(
x
)
=
F
⋅
d
i
a
g
(
x
^
)
⋅
F
H
X=C(x)=F \cdot diag(\hat{x}) \cdot F^H
X=C(x)=F⋅diag(x^)⋅FH
其中
X
X
X是循環矩陣，
x
^
\hat{x}
x^是原向量
x
x
x的傅里葉變換，
F
F
F是傅里葉變換矩陣，上標H表示共軛轉置：
X
H
=
(
X
∗
)
T
X^H=(X^{*})^T
XH=(X∗)T。<br> 換句話說，
X
X
X相似於對角陣，
X
X
X的特征值是
x
^
\hat x
x^的元素。
另一方面，如果一個矩陣能夠表示成兩個傅里葉矩陣夾一個對角陣的乘積形式，則它是一個循環矩陣。其生成向量是對角元素的傅里葉逆變換：
F
⋅
d
i
a
g
(
y
)
⋅
F
H
=
C
(
F
−
1
(
y
)
)
F \cdot diag(y) \cdot F^H = C(\mathcal{F}^{-1}(y))
F⋅diag(y)⋅FH=C(F−1(y))<br> 這個公式初看疑問很多，以下一一討論。
X
X
X是什么？
X
X
X是由原向量
x
x
x生成的循環矩陣。以矩陣尺寸
K
=
4
K=4
K=4為例。
X
=
C
(
x
)
=
[
x
1
x
2
x
3
x
4
x
4
x
1
x
2
x
3
x
3
x
4
x
1
x
2
x
2
x
3
x
4
x
1
]
X=C(x)=\begin{bmatrix} x_1 & x_2 & x_3 & x_4 \\ x_4 & x_1 & x_2 & x_3 \\ x_3 & x_4 & x_1 & x_2 \\ x_2 & x_3 & x_4 & x_1 \\ \end{bmatrix}
X=C(x)=⎣⎢⎢⎡x1x4x3x2x2x1x4x3x3x2x1x4x4x3x2x1⎦⎥⎥⎤
F
F
F 是什么？
F
F
F是離散傅里葉矩陣（DFT matrix），可以用一個復數
ω
=
e
−
2
π
i
/
K
\omega = e^{-2\pi i/K}
ω=e−2πi/K表示，其中
K
K
K為方陣
F
F
F的尺寸。以
K
=
4
K=4
K=4為例。
F
=
1
K
[
1
1
1
1
1
ω
ω
2
ω
3
1
ω
2
ω
4
ω
6
1
ω
3
ω
6
ω
9
]
F=\frac{1}{\sqrt{K}} \begin{bmatrix} 1 & 1 & 1 & 1 \\ 1 & \omega & \omega^2 & \omega^3 \\ 1 & \omega^2 & \omega ^4 & \omega^6 \\ 1 & \omega^3 & \omega^6 & \omega^9 \end{bmatrix}
F=K
<svg width="400em" height="1.08em" viewbox="0 0 400000 1080" preserveaspectratio="xMinYMin slice">
<path d="M95,702c-2.7,0,-7.17,-2.7,-13.5,-8c-5.8,-5.3,-9.5,
-10,-9.5,-14c0,-2,0.3,-3.3,1,-4c1.3,-2.7,23.83,-20.7,67.5,-54c44.2,-33.3,65.8, -50.3,66.5,-51c1.3,-1.3,3,-2,5,-2c4.7,0,8.7,3.3,12,10s173,378,173,378c0.7,0, 35.3,-71,104,-213c68.7,-142,137.5,-285,206.5,-429c69,-144,104.5,-217.7,106.5, -221c5.3,-9.3,12,-14,20,-14H400000v40H845.2724s-225.272,467,-225.272,467 s-235,486,-235,486c-2.7,4.7,-9,7,-19,7c-6,0,-10,-1,-12,-3s-194,-422,-194,-422 s-65,47,-65,47z M834 80H400000v40H845z"></path> </svg>1⎣⎢⎢⎡11111ωω2ω31ω2ω4ω61ω3ω6ω9⎦⎥⎥⎤
把
ω
\omega
ω想象成一個角度為
2
π
/
K
2\pi/K
2π/K的向量，這個矩陣的每一行是這個向量在不斷旋轉。從上到下，旋轉速度越來越快。
之所以稱為DFT matrix，是因為一個信號的DFT變換可以用此矩陣的乘積獲得：
x
^
=
D
F
T
(
x
)
=
K
⋅
F
⋅
x
\hat{x}=DFT(x)=\sqrt{K}\cdot F \cdot x
x^=DFT(x)=K
<svg width="400em" height="1.08em" viewbox="0 0 400000 1080" preserveaspectratio="xMinYMin slice">
<path d="M95,702c-2.7,0,-7.17,-2.7,-13.5,-8c-5.8,-5.3,-9.5,
-10,-9.5,-14c0,-2,0.3,-3.3,1,-4c1.3,-2.7,23.83,-20.7,67.5,-54c44.2,-33.3,65.8, -50.3,66.5,-51c1.3,-1.3,3,-2,5,-2c4.7,0,8.7,3.3,12,10s173,378,173,378c0.7,0, 35.3,-71,104,-213c68.7,-142,137.5,-285,206.5,-429c69,-144,104.5,-217.7,106.5, -221c5.3,-9.3,12,-14,20,-14H400000v40H845.2724s-225.272,467,-225.272,467 s-235,486,-235,486c-2.7,4.7,-9,7,-19,7c-6,0,-10,-1,-12,-3s-194,-422,-194,-422 s-65,47,-65,47z M834 80H400000v40H845z"></path> </svg>⋅F⋅x
反傅里葉變換也可以通過類似手段得到：
x
=
1
K
F
−
1
x
^
x=\frac{1}{\sqrt{K}}F^{-1}\hat{x}
x=K
<svg width="400em" height="1.08em" viewbox="0 0 400000 1080" preserveaspectratio="xMinYMin slice">
<path d="M95,702c-2.7,0,-7.17,-2.7,-13.5,-8c-5.8,-5.3,-9.5,
-10,-9.5,-14c0,-2,0.3,-3.3,1,-4c1.3,-2.7,23.83,-20.7,67.5,-54c44.2,-33.3,65.8, -50.3,66.5,-51c1.3,-1.3,3,-2,5,-2c4.7,0,8.7,3.3,12,10s173,378,173,378c0.7,0, 35.3,-71,104,-213c68.7,-142,137.5,-285,206.5,-429c69,-144,104.5,-217.7,106.5, -221c5.3,-9.3,12,-14,20,-14H400000v40H845.2724s-225.272,467,-225.272,467 s-235,486,-235,486c-2.7,4.7,-9,7,-19,7c-6,0,-10,-1,-12,-3s-194,-422,-194,-422 s-65,47,-65,47z M834 80H400000v40H845z"></path> </svg>1F−1x^
傅里葉矩陣有許多性質：
ω
\omega
ω的規律即可知；</li><li>滿足
F
H
F
=
F
F
H
=
I
F^HF=FF^H=I
FHF=FFH=I，也就是說它是個**酉矩陣**（unitary）。可以通過將
F
H
F^{H}
FH展開成
ω
H
\omega^{H}
ωH驗證。</li>
注意：
F
F
F是常數，可以提前計算好，和要處理的
x
x
x無關。
對角化怎么理解？
把原公式兩邊乘以逆矩陣：
F
−
1
⋅
X
⋅
(
F
H
)
−
1
=
d
i
a
g
(
x
^
)
F^{-1} \cdot X \cdot (F^H)^{-1}=diag(\hat{x})
F−1⋅X⋅(FH)−1=diag(x^)
左
邊
=
F
H
X
F
=
d
i
a
g
(
x
^
)
左邊=F^{H}XF=diag(\hat{x})
左邊=FHXF=diag(x^)
也就是說，矩陣
X
X
X通過相似變換
F
F
F變成對角陣
d
i
a
g
(
x
^
)
diag(\hat{x})
diag(x^)，即對循環矩陣
X
X
X進行對角化。<br> 另外，
F
H
X
F
F^HXF
FHXF是矩陣
X
X
X的2D DFT變換。即傅里葉變換可以把循環矩陣對角化。
怎么證明？
可以用構造特征值和特征向量的方法證明（參看這篇論文[1](#fn1)的3.1節），此處簡單描述。
X
⋅
f
k
=
x
^
k
⋅
f
k
X\cdot f_k=\hat x_k\cdot f_k
X⋅fk=x^k⋅fk
其中
f
k
f_k
fk表示DFT矩陣的第k列，
x
^
k
\hat x_k
x^k表示傅里葉變換的第k個元素。等價於求證：
f
k
f_k
fk和
x
^
k
\hat x_k
x^k是
X
X
X的一對特征向量和特征值。
左邊向量的第i個元素為：
l
e
f
t
i
=
[
x
i
,
f
k
]
left_i = \left[ x^i, f_k \right]
lefti=[xi,fk]。
x
i
x_i
xi表示把生成向量
x
x
x向右移動i位，
[
]
[]
[]表示內積。<br> 內積只和兩個向量的相對位移有關，所以可以把
f
k
f_k
fk向左移動i位：
l
e
f
t
i
=
[
x
,
f
k
−
i
]
left_i=\left[ x, f_k^{-i}\right]
lefti=[x,fk−i]。<br> DFT矩陣列的移位可以通過數乘
ω
\omega
ω的冪實現：
f
k
i
=
f
0
⋅
ω
i
k
f_k^i=f_0\cdot \omega^{ik}
fki=f0⋅ωik。
舉例：
K
=
3
K=3
K=3，<br>
F
=
1
K
[
1
1
1
1
ω
ω
2
1
ω
2
ω
4
]
F=\frac{1}{\sqrt{K}} \begin{bmatrix} 1 & 1 & 1\\ 1 & \omega & \omega^2 \\ 1 & \omega^2 & \omega ^4 \end{bmatrix}
F=K
<svg width="400em" height="1.08em" viewbox="0 0 400000 1080" preserveaspectratio="xMinYMin slice">
<path d="M95,702c-2.7,0,-7.17,-2.7,-13.5,-8c-5.8,-5.3,-9.5,
-10,-9.5,-14c0,-2,0.3,-3.3,1,-4c1.3,-2.7,23.83,-20.7,67.5,-54c44.2,-33.3,65.8, -50.3,66.5,-51c1.3,-1.3,3,-2,5,-2c4.7,0,8.7,3.3,12,10s173,378,173,378c0.7,0, 35.3,-71,104,-213c68.7,-142,137.5,-285,206.5,-429c69,-144,104.5,-217.7,106.5, -221c5.3,-9.3,12,-14,20,-14H400000v40H845.2724s-225.272,467,-225.272,467 s-235,486,-235,486c-2.7,4.7,-9,7,-19,7c-6,0,-10,-1,-12,-3s-194,-422,-194,-422 s-65,47,-65,47z M834 80H400000v40H845z"></path> </svg>1⎣⎡1111ωω21ω2ω4⎦⎤</p>
利用
ω
N
=
1
\omega^N=1
ωN=1.<br>
f
1
⋅
ω
=
[
1
,
ω
,
ω
2
]
⋅
ω
=
[
ω
,
ω
2
,
ω
3
]
=
[
ω
,
ω
2
,
1
]
=
f
1
−
1
f_1\cdot \omega =[1,\omega,\omega^2]\cdot \omega = [\omega, \omega^2,\omega^3] = [\omega, \omega^2, 1]=f_1^{-1}
f1⋅ω=[1,ω,ω2]⋅ω=[ω,ω2,ω3]=[ω,ω2,1]=f1−1</p>
f
1
⋅
ω
2
=
[
1
,
ω
,
ω
2
]
⋅
ω
2
=
[
ω
2
,
ω
3
,
ω
4
]
=
[
ω
2
,
1
,
ω
]
=
f
1
−
2
f_1\cdot \omega^2 =[1,\omega,\omega^2]\cdot \omega^2 = [\omega^2, \omega^3,\omega^4] = [\omega^2, 1, \omega]=f_1^{-2}
f1⋅ω2=[1,ω,ω2]⋅ω2=[ω2,ω3,ω4]=[ω2,1,ω]=f1−2</p>
於是有：
l
e
f
t
i
=
[
x
,
f
k
⋅
ω
i
k
]
=
[
x
,
f
k
]
⋅
ω
i
k
left_i= \left [x, f_k \cdot \omega^{ik}\right ]=\left [x, f_k\right ] \cdot \omega^{ik}
lefti=[x,fk⋅ωik]=[x,fk]⋅ωik
右邊的$\hat x = F\cdot x
，
考
慮
到
，考慮到
，考慮到F
的
第
k
行
和
第
k
列
相
同
，
的第k行和第k列相同，
的第k行和第k列相同，\hat x_k=[f_k,x]
。
另
外
。另外
。另外f_k
的
第
i
個
元
素
為
的第i個元素為
的第i個元素為\omega^{ik}$：<br>
r
i
g
h
t
i
=
x
^
k
⋅
f
k
i
=
[
f
k
∗
,
x
]
⋅
ω
k
i
right_i = \hat x_k\cdot f_{ki}=[f_k^*,x]\cdot \omega^{ki}
righti=x^k⋅fki=[fk∗,x]⋅ωki
對任意k列的第i個元素有：
l
e
f
t
i
=
r
i
g
h
t
i
left_i=right_i
lefti=righti
更多性質
利用對角化，能推導出循環矩陣的許多性質。
轉置
循環矩陣的轉置也是一個循環矩陣（可以查看循環矩陣各元素排列證明），其特征值和原特征值共軛。
X
T
=
F
⋅
d
i
a
g
(
(
x
^
)
∗
)
⋅
F
H
X^T=F \cdot diag((\hat{x})^*) \cdot F^H
XT=F⋅diag((x^)∗)⋅FH
X
T
=
(
F
H
)
T
⋅
d
i
a
g
(
x
^
)
F
T
X^T=(F^H)^T\cdot diag(\hat x) F^T
XT=(FH)T⋅diag(x^)FT<br> 由於
F
F
F是對稱酉矩陣，且已知
X
X
X是實矩陣：<br>
X
T
=
F
∗
⋅
d
i
a
g
(
x
^
)
F
=
(
F
∗
⋅
d
i
a
g
(
x
^
)
F
)
∗
=
F
⋅
d
i
a
g
(
(
x
^
)
∗
)
⋅
F
H
X^T=F^*\cdot diag(\hat x) F=\left( F^*\cdot diag(\hat x) F\right)^*=F\cdot diag((\hat x)^*) \cdot F^H
XT=F∗⋅diag(x^)F=(F∗⋅diag(x^)F)∗=F⋅diag((x^)∗)⋅FH
如果原生成向量
x
x
x是對稱向量（例如[1,2,3,4,3,2]），則其傅里葉變換為實數，則：<br>
X
T
=
C
(
F
−
1
(
x
^
)
∗
)
=
C
(
x
)
X^T=C\left( \mathcal F^{-1}(\hat x)^* \right)=C(x)
XT=C(F−1(x^)∗)=C(x)
卷積
循環矩陣乘向量等價於生成向量的逆序和該向量卷積，可進一步轉化為福利也變化相乘。
注意卷積本身即包含逆序操作，另外利用了信號與系統中經典的“時域卷積，頻域相乘”。
F
(
X
y
)
=
F
(
C
(
x
)
y
)
=
F
(
x
ˉ
∗
y
)
=
F
∗
(
x
)
⊙
F
(
y
)
\mathcal{F}(Xy) = \mathcal{F}(C(x)y)=\mathcal{F}(\bar x*y)=\mathcal{F}^*(x)\odot \mathcal{F}(y)
F(Xy)=F(C(x)y)=F(xˉ∗y)=F∗(x)⊙F(y)<br> 其中
x
ˉ
\bar x
xˉ表示把
x
x
x的元素倒序排列。星號表示共軛。
相乘
設
C
,
B
C,B
C,B為循環矩陣，其乘積的特征值等於特征值的乘積：
A
B
=
F
⋅
d
i
a
g
(
a
^
)
⋅
F
H
⋅
F
⋅
d
i
a
g
(
b
^
)
⋅
F
H
AB = F\cdot diag(\hat a) \cdot F^H \cdot F \cdot diag(\hat b) \cdot F^H
AB=F⋅diag(a^)⋅FH⋅F⋅diag(b^)⋅FH
=
F
⋅
d
i
a
g
(
a
^
)
⋅
d
i
a
g
(
b
^
)
⋅
F
H
=
F
⋅
d
i
a
g
(
a
^
⊙
b
^
)
⋅
F
H
= F\cdot diag(\hat a) \cdot diag(\hat b) \cdot F^H=F\cdot diag(\hat a \odot \hat b) \cdot F^H
=F⋅diag(a^)⋅diag(b^)⋅FH=F⋅diag(a^⊙b^)⋅FH
=
C
(
F
−
1
(
a
^
⊙
b
^
)
)
=C\left( \mathcal F^{-1}(\hat a \odot \hat b)\right)
=C(F−1(a^⊙b^))<br> 乘積也是循環矩陣，其生成向量是原生成向量對位相乘的傅里葉逆變換。
相加
和的特征值等於特征值的和：
A
+
B
=
F
⋅
d
i
a
g
(
a
^
)
⋅
F
H
+
F
⋅
d
i
a
g
(
b
^
)
⋅
F
H
=
F
⋅
d
i
a
g
(
a
^
+
b
^
)
⋅
F
H
A+B = F\cdot diag(\hat a) \cdot F^H + F \cdot diag(\hat b) \cdot F^H=F\cdot diag(\hat a +\hat b) \cdot F^H
A+B=F⋅diag(a^)⋅FH+F⋅diag(b^)⋅FH=F⋅diag(a^+b^)⋅FH
=
C
(
F
−
1
(
a
^
+
b
^
)
)
=
C
(
F
−
1
(
a
^
)
+
F
−
1
(
b
^
)
)
=
C
(
a
+
b
)
=C\left( \mathcal F^{-1} (\hat a +\hat b)\right)=C\left( \mathcal F^{-1} (\hat a )+F^{-1} (\hat b )\right)=C\left( a+b\right)
=C(F−1(a^+b^))=C(F−1(a^)+F−1(b^))=C(a+b)<br> 和也是循環矩陣，其生成向量是原生成向量的和。
求逆
循環矩陣的逆，等價於將其特征值求逆。
X
−
1
=
F
⋅
d
i
a
g
(
x
^
)
−
1
F
H
=
C
(
F
−
1
(
d
i
a
g
(
x
^
)
−
1
)
)
X^{-1}=F\cdot diag(\hat{x})^{-1}F^H=C\left( \mathcal F^{-1}(diag(\hat{x})^{-1}) \right)
X−1=F⋅diag(x^)−1FH=C(F−1(diag(x^)−1))
F
⋅
d
i
a
g
(
x
^
)
−
1
⋅
F
H
⋅
F
⋅
d
i
a
g
(
x
^
)
⋅
F
H
F\cdot diag(\hat{x})^{-1}\cdot F^H \cdot F\cdot diag(\hat{x}) \cdot F^H
F⋅diag(x^)−1⋅FH⋅F⋅diag(x^)⋅FH
=
F
⋅
d
i
a
g
(
x
^
)
−
1
⋅
d
i
a
g
(
x
^
)
⋅
F
H
=
F
⋅
F
H
=
I
=F\cdot diag(\hat{x})^{-1}\cdot diag(\hat{x}) \cdot F^H=F\cdot F^H=I
=F⋅diag(x^)−1⋅diag(x^)⋅FH=F⋅FH=I
逆也是循環矩陣
有什么用？
該性質可以將循環矩陣的許多運算轉換成更簡單的運算。例如：
X
H
X
=
F
⋅
d
i
a
g
(
x
^
⊙
x
^
∗
)
⋅
F
H
=
C
(
F
−
1
(
x
^
⊙
x
^
∗
)
)
X^HX=F \cdot diag(\hat{x} \odot \hat{x}^*)\cdot F^H=C\left( \mathcal F^{-1}(\hat{x} \odot \hat{x}^*) \right)
XHX=F⋅diag(x^⊙x^∗)⋅FH=C(F−1(x^⊙x^∗))
原始計算量：兩個方陣相乘（
O
(
K
3
)
O(K^3)
O(K3)）<br> 轉化后的計算量：反向傅里葉(
K
log
K
K\log K
KlogK)+向量點乘（
K
K
K）。
CV的許多算法中，都利用了這些性質提高運算速度，例如2015年TPAMI的這篇高速跟蹤KCF方法[2](#fn2)。
二維情況
以上探討的都是原始信號為一維的情況。以下證明二維情況下的
F
H
X
F
=
d
i
a
g
(
x
^
)
F^HXF=diag(\hat x)
FHXF=diag(x^)，推導方法和一維類似。
x
x
x是二維生成矩陣，尺寸
N
×
N
N\times N
N×N。<br>
X
X
X是一個
N
2
×
N
2
N^2\times N^2
N2×N2的分塊循環矩陣，其uv塊記為
x
u
v
x^{uv}
xuv，表示
x
x
x下移u行，右移v列。<br>
F
F
F是
N
2
×
N
2
N^2\times N^2
N2×N2的二維DFT矩陣，其第uv塊記為
f
u
v
=
{
ω
u
i
+
v
j
}
i
j
f_{uv}=\left\{ \omega^{ui+vj} \right\} _{ij}
fuv={<!-- -->ωui+vj}ij。
舉例：N=3
f
01
=
[
1
ω
ω
2
1
ω
ω
2
1
ω
ω
2
]
,
f
11
=
[
1
ω
ω
2
ω
ω
2
ω
3
ω
2
ω
3
ω
4
]
f_{01}=\begin{bmatrix} 1 & \omega & \omega^2 \\ 1 & \omega & \omega^2 \\ 1 & \omega & \omega^2\end{bmatrix}, f_{11}=\begin{bmatrix} 1 & \omega & \omega^2 \\ \omega & \omega^2 & \omega^3 \\ \omega^2 & \omega^3 & \omega^4\end{bmatrix}\\
f01=⎣⎡111ωωωω2ω2ω2⎦⎤,f11=⎣⎡1ωω2ωω2ω3ω2ω3ω4⎦⎤</p>
需要驗證的共有
N
×
N
N\times N
N×N個等式，其中第
u
v
uv
uv個為：<br>
[
X
,
f
u
v
]
=
x
^
u
v
⋅
f
u
v
[X, f_{uv}]=\hat x_{uv}\cdot f_{uv}
[X,fuv]=x^uv⋅fuv<br> 其中
[
X
,
f
u
v
]
[X, f_{uv}]
[X,fuv]表示把
x
u
v
x^{uv}
xuv分別和
f
u
v
f_{uv}
fuv做點乘，結果矩陣元素求和。
[
x
i
j
,
f
u
v
]
=
x
^
u
v
⋅
ω
u
i
+
v
j
[x^{ij},f_{uv}]=\hat x_{uv} \cdot \omega^{ui+vj}
[xij,fuv]=x^uv⋅ωui+vj
再次利用兩個性質：1) 點乘只和兩個矩陣相對位移有關，2)
f
u
v
f_{uv}
fuv的位移可以用乘
ω
\omega
ω冪實現：<br>
l
e
f
t
i
j
=
[
x
,
f
u
v
−
i
−
j
]
=
[
x
,
f
u
v
]
⋅
ω
u
i
+
v
j
=
x
^
u
v
⋅
ω
u
i
+
v
j
=
r
i
g
h
t
i
j
left_{ij}=[x,f_{uv}^{-i-j}]=[x,f_{uv}]\cdot \omega^{ui+vj}=\hat x_{uv} \cdot \omega^{ui+vj}=right_{ij}
leftij=[x,fuv−i−j]=[x,fuv]⋅ωui+vj=x^uv⋅ωui+vj=rightij
代碼
以下matlab代碼驗證上述性質。需要注意的是，matlab中的dftmatx函數給出的結果和本文定義略有不同，需做一簡單轉換。另外，matlab中的撇號表示共軛轉置，transpose為轉置函數，conj為共軛函數。
clear;clc;close all;
% 1. diagnolize
K = 5; % dimension of problem
x_base = rand(1,K); % generator vector
X = zeros(K,K); % circulant matrix
for k=1:K
X(k,:) = circshift(x_base, [0 k-1]);
end
x_hat = fft(x_base); % DFT
F = transpose(dftmtx(K))/sqrt(K); % the " ' " in matlab is transpose + conjugation
X2 = F*diag(x_hat)*F';
display(X);
display(real(X2));
% 2. fast compute correlation
C = X'*X;
C2 = (x_hat.*conj(x_hat))*conj(F)/sqrt(K);
display(C);
display(C2);