http://blog.sina.com.cn/s/blog_9206acb20101chc9.html
http://www.cnblogs.com/LeftNotEasy/archive/2011/01/19/svd-and-applications.html
http://rabbit3306.diandian.com/post/2012-06-12/40027084190
同样因为Dictionary Learning as a Non-convex Optimization Problem的珠玉在前，因此自己编出了K-SVD的字典学习方法。
K-SVD和MOD一样的地方在于，每次迭代中稀疏表示的迭代是通过追踪算法完成的。
不一样的地方在于，K-SVD在进行字典迭代的时候，不是像MOD那样用全局算法一次就算出来，而是将焦点集中在字典的每一列。每次取出字典的一列，以及X中和这一列对应的一行，算出这一列和这一行的贡献，也就是用Y-其他列*其他行，减出来的矩阵就是贡献。由于X是稀疏矩阵，在取出来的这一行中只有有限几个值是非零的，我们可以认为在贡献中只有非零值对应的列才是真的贡献。因此可以对贡献进行缩减。对得到的缩减矩阵进行SVD变换，也就是取缩减矩阵的最大近似。缩减矩阵本身应该=字典中取出的一行*X中取出的一列，因此进行SVD分解后，字典中的一行相当于u1，而稀疏表示=s(1,1)*v1；
for i=1:20 % 循环20次
X=zeros(p,m);
%%更新稀疏表示：
for j=1:m
R=Y;
xomp=zeros(p,1);
C = D'*R(:,j);
[v,PO] = sort((abs(C)),'descend');
DI=PO(1:4);
% update the coefficients
xomp(DI) = xomp(DI) + C(DI);
% perform a back projection of the coefficients
sel = find(xomp~=0);
xomp = zeros(p,1);
xomp(sel) = D(:,sel) \ Y(:,j); % Back-Project，OMP的方法
X(:,j)=xomp; % 更新了X
end
%%更新字典：
for k=1:n
[ro,omiga]= find(X(k,:)); % 找到取出的X的这一行中非零元的位置
Xt1=X(1:(k-1),:);
Xt2=X((k+1):p,:);
E = Y-D(:,1:(k-1))*Xt1-D(:,(k+1):p)*Xt2; %减去另外的行*另外的列，得到贡献矩阵
Ej=E(:,omiga); % 缩减贡献，只有对应X的那一行的非零元的列才是真的贡献
[u,s,v]=svd(Ej);
if ~isempty(s) % 有时候对应的那一行一个非零元也没有，所以要做个判断
D(:,k)=u(:,1); %更新了字典的一列
X(k,omiga)=s(1,1)*v(:,1); %更新了稀疏表示，这里的更新用来干嘛？
end
end
end
强烈推荐下这个网站http://www.cnblogs.com/LeftNotEasy/ 。比如这篇机器学习中的数学（5）SVD分解及其应用简直是太wonderful了，像《数学之美》一样有大家写小书的感觉，很赞。
Dictionary Learning as a Non-convex Optimization Problem 作为非凸优化问题的字典学习的算法解释
模拟参数的设置：
if not(exist('w'))
w = 10; % 块的边长
end
n = w*w;
p = 2*n; % 字典中原子的个数
m = 20*p; % 训练用块数的多少
k = 4; % 希望得到的稀疏度
训练样本用于产生字典：
if not(exist('f'))
f = rescale( crop(load_image('barb'),256) ); % 截取原图像的中心区域
end
n0 = size(f,1);
选取训练字典用的块数：
q = 3*m; % 多余训练用的样本，后面将会筛选
x = floor( rand(1,1,q)*(n0-w) )+1; % 随机选取块的中心像素位置，x坐标，(n0-w)是要保证块是完整的
y = floor( rand(1,1,q)*(n0-w) )+1; % 随机选取块的中心像素位置，y坐标
抽取训练用的块，并存储在Y中：
[dY,dX] = meshgrid(0:w-1,0:w-1); % 建立块中的坐标
Xp = repmat(dX,[1 1 q]) + repmat(x, [w w 1]); % 选取的块中所有点的横坐标
Yp = repmat(dY,[1 1 q]) + repmat(y, [w w 1]); % 选取的块中所有点的纵坐标
Y = f(Xp+(Yp-1)*n0); % 在256*256的一个向量中取值，坐标值都转换成向量中的形式
Y = reshape(Y, [n q]);
去除每个块中的均值：
Y = Y - repmat( mean(Y), [n 1] );
只取其中这些块中能量较大的m个（求这些块中比较白的那些，为什么？）：
[tmp,I] = sort(sum(Y.^2), 'descend'); % 按照递减排列
Y = Y(:,I(1:m)); % 取前m个
每个块是如何产生字典的：
ProjC = @(D)D ./ repmat( sqrt(sum(D.^2)), [w^2, 1] ); % 归一化
sel = randperm(m);
sel = sel(1:p);
D0 = ProjC( Y(:,sel) ); % 随机在m个块中取p个作为初始的字典（原来这里的atom指的是一个块？！）
D = D0;
定义迭代函数：
select = @(A,k)repmat(A(k,:), [size(A,1) 1]); % @为定义一个函数，相当于子函数，括号里为参数
ProjX = @(X,k)X .* (abs(X) >= select(sort(abs(X), 'descend'),k)); % 只保留X中最大的k个，保证k稀疏
迭代X：在x的0范数小于k0的情况下求解y-Dx的二次方的最小值是一个非凸非光滑的优化问题，可以用前向-后向迭代的方法求解（投影梯度减小），
for i=1:niter
R = D*X-Y;
E(end+1,:) = sum(R.^2);
X = ProjX(X - gamma * D'*R, k);
end
（y-D*x）的二次方对x求导得到梯度（D*D*x-D*y），x每次都沿着梯度走lambda步，迭代。即为投影梯度减小。
迭代D：和迭代X相同， （y-D*x）的二次方对D求导得到梯度（x*x*D-x*y），D每次都沿着梯度走tau步，迭代。
for i=1:niter_dico
R = D*X - Y;
E(end+1) = sum(R(:).^2);
D = ProjC( D + tau * (Y-D*X)*X' );
end
X和D交错迭代，最后就得到了符合k稀疏度的字典D。
5929
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
