您的位置：首页 > 其它

循环矩阵傅里叶对角化

2016-12-31 21:15 218 查看

原文：http://blog.csdn.net/shenxiaolu1984/article/details/50884830

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

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

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

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

有什么用？

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

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
[/code]

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. ↩

内容来自用户分享和网络整理，不保证内容的准确性，如有侵权内容，可联系管理员处理

标签：

相关文章推荐

新的分享

章节导航