Mar 11, 2010

matrix of linear transform and its transpose

线性变换可以用矩阵 A 表示,但是 A 可能很大,或者不好直接表示,这时在实现(Matlab)中就使用 function handle 或者 class object.

如 l1_ls 等算法中需要 A 的 Transpose (应该是 conjugate transpose?)

Ax, A^T*y
若 A=RW, 则 A^T = R^T W^T

1. A:DWT, DCT 等,参见 dftmtx, dftmtx2 等
A^T 就是反变换(或者正变换,如果A是反变换)
A^T *A = I,这要求正变换反变换的系数一致,比如很多文献中(包括matlab 中),DFT 定义为
\sum_m \sum_n f(m,n) \exp ...
IDFT 定义为
1/(MN) \sum_u \sum_v F(u,v) \exp ...
这样的定义不满足A^T A = I, 而是  A^T A =MN I
我们只需采用另一种形式即可
1/sqrt(MN) \sum_m \sum_n f(m,n) \exp ...

1/sqrt(MN) \sum_u \sum_v F(u,v) \exp...
注意:在采用function handle 或者 class object.时,还是可以采用传统定义(如 matlab fft2)的 DFT,只要记住A^T 是相应的反变换,即可。

2. A*z: subset subsampling, 取 z 的一部分,令选中的下标集为 I
z=A^T*w: z=zeros(n,1); z(I)=w

3. A: space-invariant circular (periodic) convolutions
参见论文  An EM Algorithm for Wavelet-Based Image Restoration
A 可以分解为
A=U^H D U
In the above equation, U is the matrix that represents the 2-D discrete Fourier transform, U^H is its inverse (U is an orthogonal matrix, that is, U^H U = I, where H denotes conjugate transpose), and D is a diagonal matrix containing the DFT coefficients of the convolution operator represented by A. This means that multiplication by A can be performed in the discrete Fourier domain with a simple point-wise multiplication (recall that D is diagonal).
Ax = U^H D Ux = U^H D \tilde{x}
\tilde{x}is the DFT of .

thus A^H z = U^H D^H  \tilde{z}

If matrix A is block-Toeplitz, but not block-circulant, it is possible to embed the nonperiodic convolution that it represents in a larger periodic convolution and still work in the DFT domain [16]. Accordingly, all the results and statements made in this paper concerning circulant observation matrices (periodic convolutions) can be extended to the Toeplitz case.

% % center and normalize the blur
h = fftshift(h);   
h = h/sum(h(:));

% definde the function handles that compute 
% the blur and the conjugate blur.
R = @(x) real(ifft2(fft2(h).*fft2(x)));
RT = @(x) real(ifft2(conj(fft2(h)).*fft2(x)));

参见 matlab dftmtx, convmtx (convolution matrix), convmtx2 (2D convolution matrix)
以及我写的 dftmtx2 (2D dft transform matrix)
function A = dftmtx2(M,N)
% 2D DFT transform matrix
% A = dftmtx2(M,N)
% y = Ax equal to
% Y = fft2(X), x = X'(:),y = Y'(:)
% row wise

d = M*N;
A = zeros(d,d);

[X Y]=meshgrid(0:N-1,0:M-1);
Y = Y/M;
X = X/N;
for i=0:M-1
    for j=0:N-1
        curr = exp(-2i*pi*(Y.*i + X.*j));
        curr = curr.';
        A(i*N+j+1,:)=curr(:);
    end
end
A^T * A = cI (并非等于 I)
不过A^T 确实起到了 Inverse transform 的作用

0 comments: