参考文献

非负矩阵分解(Nonnegative Matrix Factorization),简称NMF,是由Lee和Seung于1999年在自然杂志上提出的一种矩阵分解。它使分解后的所有分量均为非负值(要求纯加性的描述)并且同时实现非线性的维数约减。

非负矩阵分解将一个矩阵 $V$ 近似分解成两个非负的矩阵 $W$ 和 $H$，即$V \approx W \times H$，其中$W,H \geq 0$。那么利用迭代优化的思想使得 $W \times H$ 的结果尽可能接近$V$，有

$$ \min_{W,H} f(W,H) = \frac{1}{2} \sum_{i=1}^n \sum_{j=1}^m (V_{ij} - (W \times H)_{ij})^2 \\ \mathsf{s.t.}\ W_{ia} \geq 0, H_{bj} \geq 0, \forall i,a,b,j $$

其中$s.t.$（subject to）为“受限于，服从于”的意思。

标准方法

$$ W^{k+1} = \max(0,W^k - \alpha \cdot \nabla_Wf(W^k,H^k)) \\ H^{k+1} = \max(0,H^k - \alpha \cdot \nabla_Hf(W^k,H^k)) $$

上述方法遵从标准化有界优化问题的求解，即对于标准化问题

$$ \min_{x \in R^n} f(x) \\ \mathsf{s.t.}\ l_i \leq x_i \leq u_i, i = 1, 2, ...,n $$

其解为

$$ x^{k+1} = P[(x^k - \alpha \cdot \nabla f(x^k))] $$

其中

$$ P[x_i] = \begin{cases} x_i & \mathsf{if}\ l_i < x_i < u_i, \\ l_i & \mathsf{if}\ x_i \leq l_i, \\ u_i & \mathsf{if}\ x_i \geq u_i, \end{cases} $$

乘法更新法

1999年，Lee和Seung在自然杂志上提出的一种矩阵分解方法，该方法收敛速更快。

$$ W_{ia}^{k+1} = W_{ia}^{k} \frac{(V(H^k)^\mathsf{T})_{ia}}{(W^kH^k(H^k)^\mathsf{T})_{ia}} \ \forall i,a \\ H_{bj}^{k+1} = H_{bj}^{k} \frac{((W^{k+1})^\mathsf{T}V)_{bj}}{((W^{k+1})^\mathsf{T} W^{k+1} H^k)_{bj}} \ \forall b,j $$

还有进一步改进的方法

详见《Projected Gradient Methods for Nonnegative Matrix Factorization》

示例

以下展示了非负矩阵分解在图像领域的应用示例。

V=double(imread('图片名或路径'));      #	这里的图片默认为RGB格式（即三通道图片）
imshow(mat2gray(V));
V1=V(:,:,1);		# 将图片的第一通道数值赋值给V1
V2=V(:,:,2);		# 将图片的第二通道数值赋值给V2
V3=V(:,:,3);		# 将图片的第三通道数值赋值给V3

[i,u]=size(V1);		# 将V1的行数赋值给i，列数赋值给u
r=100;				# 设置分解矩阵的秩
W=rand(i,r);		# 初始化WH，为非负数
H=rand(r,u);		# 初始化WH，为非负数
maviter=600;		# 最大迭代次数
-----------------------------------------------------------------------------------
# 以下为非负矩阵算法公式，V2，V3同上处理，就不重复展示,以'...'代替
for iter=1:maviter
    W=W.*((V1./(W*H))*H');           	
    W=W./(ones(i,1)*sum(W));    
    H=H.*(W'*(V1./(W*H)));
end
V1=W*H;
...
-----------------------------------------------------------------------------------
# 以下程序用于重新拼接图形
img_V(:,:,1)=V1(:,:);
img_V(:,:,2)=V2(:,:);
img_V(:,:,3)=V3(:,:);
figure;
imshow(mat2gray(img_V));

以下为改进后的程序，增加了图像的后期处理

%------------加载数据
picture = imread('airplane.png');

V=double(picture); %将图像转换成像素矩阵

%绘制原始图像
figure;
imshow(mat2gray(V));
title("原图");

H = fspecial('unsharp');
sharpened = imfilter(mat2gray(V),H,'replicate');
figure;
imshow(sharpened);
title('对原图进行锐化处理')

%------------进行非负矩阵分解
V1=V(:,:,1);
V2=V(:,:,2);
V3=V(:,:,3);

[W1,H1] = NMF(V1);
V1_NMF=W1*H1;

[W2,H2] = NMF(V2);
V2_NMF=W2*H2;

[W3,H3] = NMF(V3);
V3_NMF=W3*H3;

img_V(:,:,1)=V1_NMF(:,:);
img_V(:,:,2)=V2_NMF(:,:);
img_V(:,:,3)=V3_NMF(:,:);

%绘制分解后的图像效果
figure;
imshow(mat2gray(img_V));
title("NMF");

H = fspecial('unsharp');
sharpened = imfilter(mat2gray(img_V),H,'replicate');
figure;
imshow(sharpened);
title('对NMF图像进行锐化处理')

%------------对分解后的图像进行加强
RGB=mat2gray(img_V);
HSV = rgb2hsv(RGB);
%H=HSV(:,:,1);
%S=HSV(:,:,2);
V=HSV(:,:,3);

%imadjust函数处理后的图像
V= imadjust(V);
HSV(:,:,3)=V;
RGB_1 =  hsv2rgb(HSV);
figure;
imshow(RGB_1);
title('imadjust函数处理后的图像')

H = fspecial('unsharp');
sharpened = imfilter(RGB_1,H,'replicate');
figure;
imshow(sharpened);
title('对imadjust函数处理后的图像进行锐化处理')

%adapthisteq函数处理后的图像
V= adapthisteq(V);
HSV(:,:,3)=V;
RGB_2 =  hsv2rgb(HSV);
figure;
imshow(RGB_2);
title('adapthisteq函数处理后的图像')

H = fspecial('unsharp');
sharpened = imfilter(RGB_2,H,'replicate');
figure;
imshow(sharpened);
title('对adapthisteq函数处理后的图像进行锐化处理')

%histeq函数处理后的图像
V= histeq(V);
HSV(:,:,3)=V;
RGB_3 =  hsv2rgb(HSV);
figure;
imshow(RGB_3);
title('histeq函数处理后的图像')

H = fspecial('unsharp');
sharpened = imfilter(RGB_3,H,'replicate');
figure;
imshow(sharpened);
title('对histeq函数处理后的图像进行锐化处理')

%{
J = mat2gray(img_V);
I=rgb2gray(J);
%J1 = medfilt2(I);
img_V(:,:,3) = adapthisteq(img_V(:,:,3));
figure;
subplot(1,2,1);
imshow(I);
subplot(1,2,2);
imshow(mat2gray(img_V));
%}

%------------------------
%{
img = J;
figure;
imshow (img) ;
title("原图");
%转hsv
img_hsv = rgb2hsv(img);
%提取亮度
img_v = img_hsv(:,:,3);

%先进行直方图均衡化处理
img_hist = histeq(img_v);
%qqq = img_hist(:,:,3);

%再进行锐化处理
%w = fspecial('log',[5,5], 0.4);%生成高斯-拉普拉斯滤波器
w = fspecial('unsharp');
img_ruihua = imfilter(img_hist,w,'replicate') ;
%img_ruihua = filter(img_hist,w,5) ;
%img_ruihua = medfilt2(img_hist(:,:,3),[5,5]);

img_hsv(:,:,3) = img_ruihua;
img_after = hsv2rgb(img_hsv);
imshow(img_after);
title( '处理后')

%----------------
%锐化函数
function img_result = filter(img_v,w,muban_size)
moban_size = muban_size;
[M,N] = size(img_v);
img_lap = zeros(M,N);
expand_img = double(wextend('2D','zpd',img_v,floor(moban_size/2))); %括号扩展0，转double为了矩阵运算

for i=1:M
    for j=1:N
        ave = sum( sum( expand_img (i:i+moban_size-1,j:j+moban_size-1).* w));%号取出扩展元素与模板相乘，并求
        img_lap(i,j)= ave;
    end
end

img_result = img_v - img_lap;
end

%%直方图均衡化函数
function hist_img = hist_1(I)
I = floor(I.*255);
[M, N] = size(I);
size_img = M*N;
c=zeros(1,256);%统计每个每个像素值的个数
b= c;%转化前后的对照表
temp = I(:) ;
temp = sort (temp) ;

for i = 1:size_img
    c(temp(i)+1) = c(temp(1)+1)+1;
end
a=c;%木和
for i = 2:256
    a(i) = c(i) + a(i-1);
end

min_cdf = 10000;
for i = 1:256
    if a(i)>0
        if a(i) < min_cdf
            min_cdf= a(i);
        end
    end
end
for j = 1:256
    if a(j) > 0
        b(j) = round (255* (a(j) -min_cdf)/(size_img-min_cdf));
    end
end
for i = 1:M
    for j=1:N
        I(i,j) = b(I(i,j)+1) ;
    end
end
hist_img = I./255;
end
%}
function [W,H] = NMF(V)
[i,u]=size(V);                                    %计算V的规格
r=100;                                  %设置分解矩阵的秩
W=rand(i,r);                            %初始化WH，为非负数
H=rand(r,u);
maviter=600;                                    %最大迭代次数
for iter=1:maviter
    W=W.*((V./(W*H))*H');           %注意这里的三个公式和文中的是对应的
    W=W./(ones(i,1)*sum(W));
    H=H.*(W'*(V./(W*H)));
end
end

非负矩阵分解 / Nonnegative Matrix Factorization (NMF)

标准方法

乘法更新法

还有进一步改进的方法

示例