循环矩阵傅里叶对角化

最新推荐文章于 2024-05-14 22:15:00 发布

shenxiaolu1984

最新推荐文章于 2024-05-14 22:15:00 发布

阅读量8w

点赞数 140

分类专栏：数学文章标签：循环矩阵傅里叶变换对角化快速计算

本文链接：https://blog.csdn.net/shenxiaolu1984/article/details/50884830

版权

数学专栏收录该内容

11 篇文章 20 订阅

订阅专栏

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

文献中，一般用如下方式表达这一概念：
$\cdot diag(\hat{x}) \cdot F^H$

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

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

$\cdot diag(y) \cdot F^H = C(\mathcal{F}^{-1}(y))$
这个公式初看疑问很多，以下一一讨论。

$X$ 是什么？

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

$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}$

$F$ 是什么？

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

$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}$

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

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

$\hat{x}=DFT(x)=\sqrt{K}\cdot F \cdot x$

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

$x=\frac{1}{\sqrt{K}}F^{-1}\hat{x}$

傅里叶矩阵有许多性质：

是对称矩阵，观察 $\omega$ 的规律即可知；
满足 $F^HF=FF^H=I$ ，也就是说它是个酉矩阵（unitary）。可以通过将 $F^{H}$ 展开成 $\omega^{H}$ 验证。

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

对角化怎么理解？

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

$F^{-1} \cdot X \cdot (F^H)^{-1}=diag(\hat{x})$
利用前述酉矩阵性质：
$左边=F^{H}XF=diag(\hat{x})$

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

怎么证明？

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

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

左边向量的第i个元素为： $left_i = \left[ x^i, f_k \right]$ 。 $x_i$ 表示把生成向量 $x$ 向右移动i位， $[]$ 表示内积。
内积只和两个向量的相对位移有关，所以可以把 $f_k$ 向左移动i位： $left_i=\left[ x, f_k^{-i}\right]$ 。
DFT矩阵列的移位可以通过数乘 $\omega$ 的幂实现： $f_k^i=f_0\cdot \omega^{ik}$ 。

举例： $K = 3$ ，
$F=\frac{1}{\sqrt{K}} \begin{bmatrix} 1 & 1 & 1\\ 1 & \omega & \omega^2 \\ 1 & \omega^2 & \omega ^4 \end{bmatrix}$

利用 $\omega^N=1$ .
$f_1\cdot \omega =[1,\omega,\omega^2]\cdot \omega = [\omega, \omega^2,\omega^3] = [\omega, \omega^2, 1]=f_1^{-1}$

$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}$

于是有：
$left_i= \left [x, f_k \cdot \omega^{ik}\right ]=\left [x, f_k\right ] \cdot \omega^{ik}$

右边的$\hat x = F\cdot x $，考虑到$ F $的第 k 行和第 k 列相同，$ \hat x_k=[f_k,x] $。另外$ f_k $的第 i 个元素为$ \omega^{ik}$：
$right_i = \hat x_k\cdot f_{ki}=[f_k^*,x]\cdot \omega^{ki}$

对任意k列的第i个元素有： $left_i=right_i$

有什么用？

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

$X^HX=F \cdot diag(\hat{x} \odot \hat{x}^*)\cdot F^H=C\left( \mathcal F^{-1}(\hat{x} \odot \hat{x}^*) \right)$

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

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

二维情况

以上探讨的都是原始信号为一维的情况。以下证明二维情况下的 $F^HXF=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_{uv}=\left\{ \omega^{ui+vj} \right\} _{ij}$ 。

举例：N=3 $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}\\$

需要验证的共有 $N\times N$ 个等式，其中第 $u v$ 个为：
$f_{uv}]=\hat x_{uv}\cdot f_{uv}$
其中 $X, f_{uv}]$ 表示把 $x^{uv}$ 分别和 $f_{uv}$ 做点乘，结果矩阵元素求和。
这个等式的第ij元素为：
$[x^{ij},f_{uv}]=\hat x_{uv} \cdot \omega^{ui+vj}$

再次利用两个性质：1) 点乘只和两个矩阵相对位移有关，2) $f_{uv}$ 的位移可以用乘 $\omega$ 幂实现：
$left_{ij}=[x,f_{uv}^{-i-j}]=[x,f_{uv}]\cdot \omega^{ui+vj}=\hat x_{uv} \cdot \omega^{ui+vj}=right_{ij}$

代码

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

Gray, Robert M. Toeplitz and circulant matrices: A review. now publishers inc, 2006. ↩︎
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. ↩︎