循环矩阵傅里叶对角化

参考：http://blog.youkuaiyun.com/shenxiaolu1984/article/details/50884830

All circulant matrices are made diagonal by the Discrete Fourier Transform (DFT), regardless of the generating vector x.

任意循环矩阵可以被傅里叶变换矩阵对角化。

diag就是对角化。

文献中，一般用如下方式表达这一概念：

X=C(x)=F⋅diag(x^)⋅FH

其中 X 是循环矩阵， x^ 是原向量 x 的傅里叶变换， F 是傅里叶变换矩阵，上标H表示共轭转置： XH=(X∗)T 。
换句话说， X 相似于对角阵， X 的特征值是 x^ 的元素。

另一方面，如果一个矩阵能够表示成两个傅里叶矩阵加一个对角阵的乘积形式，则它是一个循环矩阵。其生成向量是对角元素的傅里叶逆变换：

F⋅diag(y)⋅FH=C(F−1(y))

这个公式初看疑问很多，以下一一讨论。

X 是什么？

X 是由原向量 x 生成的循环矩阵。以矩阵尺寸 K=4 为例。

X=C(x)=⎡⎣⎢⎢⎢x1x4x3x2x2x1x4x3x3x2x1x4x4x3x2x1⎤⎦⎥⎥⎥

F 是什么？

F 是离散傅里叶矩阵（DFT matrix），可以用一个复数 ω=e−2πi/K 表示，其中 K 为方阵 F 的尺寸。以 K=4 为例。

F=1K−−√⎡⎣⎢⎢⎢⎢11111ωω2ω31ω2ω4ω61ω3ω6ω9⎤⎦⎥⎥⎥⎥

把 ω 想象成一个角度为 2π/K 的向量，这个矩阵的每一行是这个向量在不断旋转。从上到下，旋转速度越来越快。

之所以称为DFT matrix，是因为一个信号的DFT变换可以用此矩阵的乘积获得：

x^=DFT(x)=K−−√⋅F⋅x

反傅里叶变换也可以通过类似手段得到：

x=1K−−√F−1x^

傅里叶矩阵有许多性质：
- 是对称矩阵，观察 ω 的规律即可知；
- 满足 FHF=FFH=I ，也就是说它是个酉矩阵（unitary）。可以通过将 FH 展开成 ωH 验证。

注意： F 是常数，可以提前计算好，和要处理的 x 无关。

对角化怎么理解？

把原公式两边乘以逆矩阵：

F−1⋅X⋅(FH)−1=diag(x^)

利用前述酉矩阵性质：

左边=FHXF=diag(x^)

也就是说，矩阵 X 通过相似变换 F 变成对角阵 diag(x^) ，即对循环矩阵 X 进行对角化。
另外， FHXF 是矩阵 X 的2D DFT变换。即傅里叶变换可以把循环矩阵对角化。

怎么证明？

可以用构造特征值和特征向量的方法证明（参看这篇论文1的3.1节），此处简单描述。
考察待证明等式的第k列：

X⋅fk=x^k⋅fk

其中 fk 表示DFT矩阵的第k列， x^k 表示傅里叶变换的第k个元素。等价于求证： fk 和 x^k 是 X 的一对特征向量和特征值。

左边向量的第i个元素为： lefti=[xi,fk] 。 xi 表示把生成向量 x 向右移动i位， [] 表示内积。
内积只和两个向量的相对位移有关，所以可以把 fk 向左移动i位： lefti=[x,f−ik] 。
DFT矩阵列的移位可以通过数乘 ω 的幂实现： fik=f0⋅ωik 。

举例： K=3 ，

F=1K−−√⎡⎣⎢1111ωω21ω2ω4⎤⎦⎥

利用 ωN=1 .

f1⋅ω=[1,ω,ω2]⋅ω=[ω,ω2,ω3]=[ω,ω2,1]=f−11

f1⋅ω2=[1,ω,ω2]⋅ω2=[ω2,ω3,ω4]=[ω2,1,ω]=f−21

于是有：

lefti=[x,fk⋅ωik]=[x,fk]⋅ωik

右边的 x^=F⋅x ，考虑到 F 的第k行和第k列相同， x^k=[fk,x] 。另外 fk 的第i个元素为 ωik ：

righti=x^k⋅fki=[f∗k,x]⋅ωki

对任意k列的第i个元素有： lefti=righti

更多性质

利用对角化，能推导出循环矩阵的许多性质。

转置

循环矩阵的转置也是一个循环矩阵（可以查看循环矩阵各元素排列证明），其特征值和原特征值共轭。

XT=F⋅diag((x^)∗)⋅FH

可以通过如下方式证明：

XT=(FH)T⋅diag(x^)FT

由于 F 是对称酉矩阵，且已知 X 是实矩阵：

XT=F∗⋅diag(x^)F=(F∗⋅diag(x^)F)∗=F⋅diag((x^)∗)⋅FH

如果原生成向量 x 是对称向量（例如[1,2,3,4,3,2]），则其傅里叶变换为实数，则：

XT=C(F−1(x^)∗)=C(x)

卷积

循环矩阵乘向量等价于生成向量的逆序和该向量卷积，可进一步转化为福利也变化相乘。
注意卷积本身即包含逆序操作，另外利用了信号与系统中经典的“时域卷积，频域相乘”。

F(Xy)=F(C(x)y)=F(x¯∗y)=F∗(x)⊙F(y)

其中 x¯ 表示把 x 的元素倒序排列。星号表示共轭。

相乘

设 C,B 为循环矩阵，其乘积的特征值等于特征值的乘积：

AB=F⋅diag(a^)⋅FH⋅F⋅diag(b^)⋅FH

=F⋅diag(a^)⋅diag(b^)⋅FH=F⋅diag(a^⊙b^)⋅FH

=C(F−1(a^⊙b^))

乘积也是循环矩阵，其生成向量是原生成向量对位相乘的傅里叶逆变换。

相加

和的特征值等于特征值的和：

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)

和也是循环矩阵，其生成向量是原生成向量的和。

求逆

循环矩阵的逆，等价于将其特征值求逆。

X−1=F⋅diag(x^)−1FH=C(F−1(diag(x^)−1))

对角阵求逆等价于对角元素求逆。以下证明：

F⋅diag(x^)−1⋅FH⋅F⋅diag(x^)⋅FH

=F⋅diag(x^)−1⋅diag(x^)⋅FH=F⋅FH=I

逆也是循环矩阵

有什么用？

该性质可以将循环矩阵的许多运算转换成更简单的运算。例如：

XHX=F⋅diag(x^⊙x^∗)⋅FH=C(F−1(x^⊙x^∗))

原始计算量：两个方阵相乘（ O(K3) ）
转化后的计算量：反向傅里叶( KlogK )+向量点乘（ K ）。

CV的许多算法中，都利用了这些性质提高运算速度，例如2015年TPAMI的这篇高速跟踪KCF方法2。

二维情况

以上探讨的都是原始信号为一维的情况。以下证明二维情况下的 FHXF=diag(x^) ，推导方法和一维类似。

x 是二维生成矩阵，尺寸 N×N 。
X 是一个 N2×N2 的分块循环矩阵，其uv块记为 xuv ，表示 x 下移u行，右移v列。
F 是 N2×N2 的二维DFT矩阵，其第uv块记为 fuv={ωui+vj}ij 。

举例：N=3

f01=⎡⎣⎢111ωωωω2ω2ω2⎤⎦⎥,f11=⎡⎣⎢1ωω2ωω2ω3ω2ω3ω4⎤⎦⎥

需要验证的共有 N×N 个等式，其中第 uv 个为：

[X,fuv]=x^uv⋅fuv

其中 [X,fuv] 表示把 xuv 分别和 fuv 做点乘，结果矩阵元素求和。
这个等式的第ij元素为：

[xij,fuv]=x^uv⋅ωui+vj

再次利用两个性质：1) 点乘只和两个矩阵相对位移有关，2) fuv 的位移可以用乘 ω 幂实现：

leftij=[x,f−i−juv]=[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);
  
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
  
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

---

1. Gray, Robert M. Toeplitz and circulant matrices: A review. now publishers inc, 2006.↩
2. Henriques, João F., et al. “High-speed tracking with kernelized correlation filters.” Pattern Analysis and Machine Intelligence, IEEE Transactions on 37.3 (2015): 583-596.↩