循环矩阵傅里叶对角化

转载至：http://blog.csdn.net/shenxiaolu1984/article/details/50884830

循环矩阵傅里叶对角化

标签： 循环矩阵傅里叶变换对角化快速计算

2016-03-15 14:10 778人阅读 评论(9) 收藏 举报

本文章已收录于：

版权声明：本文为博主原创文章，未经博主允许不得转载。

目录(?)[+]

1. XX是什么
2. FF 是什么
3. 对角化怎么理解
4. 怎么证明
5. 更多性质
   1. 转置
   2. 卷积
   3. 相乘
   4. 相加
   5. 求逆
6. 有什么用
7. 二维情况
8. 代码

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

更多性质

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

转置

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

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为共轭函数。

<code class="language-matlab hljs  has-numbering">clear;clc;close all;

<span class="hljs-comment">% 1. diagnolize </span>
K = <span class="hljs-number">5</span>;      <span class="hljs-comment">% dimension of problem</span>

x_base = <span class="hljs-built_in">rand</span>(<span class="hljs-number">1</span>,K);     <span class="hljs-comment">% generator vector</span>
X = <span class="hljs-built_in">zeros</span>(K,K);         <span class="hljs-comment">% circulant matrix</span>
<span class="hljs-keyword">for</span> k=<span class="hljs-number">1</span>:K
    X(k,:) = <span class="hljs-built_in">circshift</span>(x_base, <span class="hljs-matrix">[<span class="hljs-number">0</span> k-<span class="hljs-number">1</span>]</span>);
<span class="hljs-keyword">end</span>

x_hat = fft(x_base);    <span class="hljs-comment">% DFT</span>

F = transpose(dftmtx(K))/<span class="hljs-built_in">sqrt</span>(K);       <span class="hljs-comment">% the " ' " in matlab is transpose + conjugation</span>

X2 = F*<span class="hljs-built_in">diag</span>(x_hat)*<span class="hljs-transposed_variable">F'</span>;

display(X);
display(<span class="hljs-built_in">real</span>(X2));

<span class="hljs-comment">% 2. fast compute correlation</span>
C = <span class="hljs-transposed_variable">X'</span>*X;
C2 = (<span class="hljs-transposed_variable">x_hat.</span>*<span class="hljs-built_in">conj</span>(x_hat))*<span class="hljs-built_in">conj</span>(F)/<span class="hljs-built_in">sqrt</span>(K);

display(C);
display(C2);</code><ul class="pre-numbering" style=""><li>1</li><li>2</li><li>3</li><li>4</li><li>5</li><li>6</li><li>7</li><li>8</li><li>9</li><li>10</li><li>11</li><li>12</li><li>13</li><li>14</li><li>15</li><li>16</li><li>17</li><li>18</li><li>19</li><li>20</li><li>21</li><li>22</li><li>23</li><li>24</li><li>25</li><li>26</li></ul><ul class="pre-numbering" style=""><li>1</li><li>2</li><li>3</li><li>4</li><li>5</li><li>6</li><li>7</li><li>8</li><li>9</li><li>10</li><li>11</li><li>12</li><li>13</li><li>14</li><li>15</li><li>16</li><li>17</li><li>18</li><li>19</li><li>20</li><li>21</li><li>22</li><li>23</li><li>24</li><li>25</li><li>26</li></ul>

---

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

顶

0

踩

0

 

 

• 上一篇核方法(kernel method)超简说明
• 下一篇KCF高速跟踪详解