循环矩阵傅里叶对角化

转载自https://blog.csdn.net/shenxiaolu1984/article/details/50884830

相关KCF跟踪算法讲解：

https://blog.csdn.net/shenxiaolu1984/article/details/50905283
https://www.cnblogs.com/YiXiaoZhou/p/5925019.html

Matlab代码实现：
https://blog.csdn.net/bitopyx/article/details/81948875?spm=1001.2014.3001.5502

https://blog.csdn.net/weixin_38128100/article/details/95729653

All circulant matrices are made diagonal by the Discrete Fourier Transform (DFT), regardless of the generating vector x.
任意循环矩阵可以被傅里叶变换矩阵对角化。

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

其中$X$是循环矩阵，$\hat{x}$是原向量$x$的傅里叶变换，$F$是傅里叶变换矩阵，上标H表示共轭转置：$X^H=(X^{\ast }{)}^T$。

换句话说，$X$相似于对角阵，$X$的特征值是$\hat{x}$的元素。

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

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

$X$是什么？

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

$F$是什么？

$F$是离散傅里叶矩阵（DFT matrix），可以用一个复数$w=e^{-2\pi i/K}$表示，其中$K$为方阵$F$的尺寸，以$K=4$为例。

把$w$想象成一个角度为$2\pi /K$的向量，这个矩阵的每一行是这个向量不断旋转。从上到下，旋转速度越来越快。
之所以称为DFT matrix，是因为一个信号的DFT变换可以用此矩阵的乘积获得：

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

傅里叶矩阵有许多性质：

• 是对称矩阵，观察$w$的规律即可知；
• 满足$F^{H}F=F{F}^H=I$，也就是说它是个酉矩阵。
  注意：$F$是常数，可以提前计算好，和要处理的$x$无关。

对角化怎么理解？

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

利用酉矩阵性质：

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

怎么证明？

可以用构造特征值和特征向量的方法证明（参看这篇论文（Gray, Robert M. Toeplitz and circulant matrices: A review. now publishers inc, 2006.）的3.1节）。此处简单描述。
考察待证明等式的第k列：

其中，$f_k$表示DFT矩阵的第k列，${\hat{x}}_k$表示傅里叶变换的第k个元素。等价于求证：$f_k$和${\hat{x}}_k$是$X$的一对特征向量和特征值。
左边向量的第i个元素为：$lef{t}_i=[x^i,f_k]$。$x^i$表示把生成向量$x$向右移动i位，[]表示内积。内积只和两个向量的相对位移有关，所以可以把$f_k$向左移动i位：$lef{t}_i=[x,f_k^{-i}]$。DFT矩阵列的移位可以通过数乘$w$的幂实现：$f_k^i=f_0\cdot w^{ik}$。

更多性质

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

转置

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

可以通过如下方式证明：

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

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

卷积

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

其中$\bar{x}$表示把$x$的元素倒序排列。星号表示共轭。

相乘

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

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

相加

和的特征值等于特征值的和。

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

求逆

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

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

逆也是循环矩阵。

有什么用？

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

原始计算量：两个方阵相乘$(O(K^3))$
转化后的计算量：反向傅里叶$(KlogK)+$向量点乘$(K)$。
CV的许多算法中，都利用了这些性质提高运算速度，例如2015年TPAMI的这篇高速跟踪KCF方法。

二维情况

以上探讨的都是原始信号为一维的情况，以下证明二维情况下的$F^{H}XF=diag(\hat{x})$，推到方法和一维类似。
$x$是二维生成矩阵，尺寸$N\times N$。
$X$是一个$N^2\times N^2$的分块循环矩阵，其uv块记为$x^{uv}$，表示$x$下移u行，右移v列。
$F$是$N^2\times N^2$的二维DFT矩阵，其第uv块记为 f u v = { w u i + v j } i j f_{uv}=\{w^{ui+vj}\}_{ij} fuv​={wui+vj}ij​。

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

其中，$[X,f_{uv}]$表示把$x^{uv}$分别和$f_{uv}$做点乘，结果矩阵元素求和。
这个等式的第ij元素为：

再次利用两个性质：（1）点乘只和两个矩阵相对位移有关（2）$f_{uv}$的位移可以用乘$w$幂实现：

代码

以下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);