循环矩阵傅里叶对角化
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);