Orthogonal Least Squares （OLS）演算法流程

實驗

要利用python實現，電腦必須安裝以下程式

• python （本文用的python版本為3.5.1）
• numpy python包（本文用的版本為1.10.4）
• scipy python包（本文用的版本為0.17.0）
• pillow python包（本文用的版本為3.1.1）

python程式碼

#coding: utf-8
'''
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# DCT基作為稀疏基，重建演算法為OMP演算法 ，影象按列進行處理
# email:
 
[email protected], 
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
'''
#匯入整合庫
import math

# 匯入所需的第三方庫檔案
import  numpy as np    #對應numpy包
from PIL import Image  #對應pillow包


#讀取影象，並變成numpy型別的 array
im = np.array(Image.open('lena.bmp'))
#print (im.shape, im.dtype)uint8

#生成高斯隨機測量矩陣
sampleRate=0.7
 
  #取樣率
Phi=np.random.randn(256*sampleRate,256)

#生成稀疏基DCT矩陣
mat_dct_1d=np.zeros((256,256))
v=range(256)
for k in range(0,256):  
    dct_1d=np.cos(np.dot(v,k*math.pi/256))
    if k>0:
        dct_1d=dct_1d-np.mean(dct_1d)
    mat_dct_1d[:,k]=dct_1d/np.linalg.norm(dct_1d)

#隨機測量
img_cs_1d=np.dot(Phi,im)

#OLS演算法函式（未完成待修改）
def cs_ols(y,D):
 
    
    L=math.floor(3*(y.shape[0])/4)
    residual=y  #初始化殘差
    index=np.zeros((L),dtype=int)
    for i in range(L):
        index[i]= -1
    result=np.zeros((256))
    for j in range(L):  #迭代次數
        pos_temp=np.argsort(np.fabs(np.dot(np.linalg.pinv(D[:,pos]),y)))#對應步驟2

        product=np.fabs(np.dot(D.T,residual))
        pos=np.argmax(product)  #最大投影係數對應的位置        
        index[j]=pos
        my=np.linalg.pinv(D[:,index>=0])        
        a=np.dot(my,y)      
        residual=y-np.dot(D[:,index>=0],a)
    result[index>=0]=a
    return  result

#重建
sparse_rec_1d=np.zeros((256,256))   # 初始化稀疏係數矩陣    
Theta_1d=np.dot(Phi,mat_dct_1d)   #測量矩陣乘上基矩陣
for i in range(256):
    print('正在重建第',i,'列。。。')
    column_rec=cs_ols(img_cs_1d[:,i],Theta_1d)  #利用OMP演算法計算稀疏係數
    sparse_rec_1d[:,i]=column_rec;        
img_rec=np.dot(mat_dct_1d,sparse_rec_1d)          #稀疏係數乘上基矩陣

#顯示重建後的圖片
image2=Image.fromarray(img_rec)
image2.show()

matlab程式碼

function Demo_CS_OLS()
%------------ read in the image --------------
img=imread('lena.bmp');     % testing image
img=double(img);
[height,width]=size(img);
%------------ form the measurement matrix and base matrix ---------------
Phi=randn(floor(0.5*height),width);  % only keep one third of the original data  
Phi = Phi./repmat(sqrt(sum(Phi.^2,1)),[floor(0.5*height),1]); % normalize each column

mat_dct_1d=zeros(256,256);  % building the DCT basis (corresponding to each column)
for k=0:1:255 
    dct_1d=cos([0:1:255]'*k*pi/256);
    if k>0
        dct_1d=dct_1d-mean(dct_1d); 
    end;
    mat_dct_1d(:,k+1)=dct_1d/norm(dct_1d);
end

%--------- projection ---------
img_cs_1d=Phi*img;          % treat each column as a independent signal

%-------- recover using omp ------------
sparse_rec_1d=zeros(height,width); % height*width的0矩陣    
Theta_1d=Phi*mat_dct_1d;%測量矩陣乘上基矩陣
for i=1:width
    column_rec=OLS(img_cs_1d(:,i),Theta_1d);%演算法的目的是得到稀疏係數
    sparse_rec_1d(:,i)=column_rec';           % sparse representation 稀疏係數
end
img_rec_1d=mat_dct_1d*sparse_rec_1d;          % inverse transform 稀疏係數乘上基矩陣


%------------ show the results --------------------
figure(1)
subplot(2,2,1),imshow(uint8(img)),title('original image')
subplot(2,2,2),imagesc(Phi),title('measurement mat')
subplot(2,2,3),imagesc(mat_dct_1d),title('1d dct mat')
psnr = 20*log10(255/sqrt(mean((img(:)-img_rec_1d(:)).^2)));
subplot(2,2,4),imshow(uint8(img_rec_1d));
title(strcat('PSNR=',num2str(psnr),'dB'));

%************************************************************************%
function [s, residual] = OLS(y, A, err)

% Orthogonal Least Squares [1] for Sparse Signal Reconstruction

% Input
% A = N X d dimensional measurment matrix
% y = N dimensional observation vector
% m = sparsity of the underlying signal

% Output
% s = estimated sparse signal
% r = residual 

% [1] T. Blumensath, M. E. Davies; "On the Difference Between Orthogonal 
% Matching Pursuit and Orthogonal Least Squares", manuscript 2007

if nargin < 3
     err = 1e-5;
end

n1=length(y);
m=floor(3*n1/4);  

s = zeros(size(A,2),1);
r(:,1) = y; L = []; Psi = [];
normA=(sum(A.^2,1)).^0.5;
NI = 1:size(A,2);

for i = 2:m+1
    DR = A'*r(:,i-1);
    [v, I] = max(abs(DR(NI))./normA(NI)');
    INI = NI(I);
    L = [L' INI']';
    NI(I)=[];
    Psi = A(:,L);
    x = Psi\y;
    yApprox = Psi*x;
    r(:,i) = y - yApprox;
    if norm(r(:,end)) < err
        break;
    end
end
s(L) = x;
residual = r(:,end);

參考文章

1、Blumensath T, Davies M E. On the difference between orthogonal matching pursuit and orthogonal least squares[J]. 2007.
2、Hashemi A, Vikalo H. Sparse Linear Regression via Generalized Orthogonal Least-Squares[J]. arXiv preprint arXiv:1602.06916, 2016.

歡迎python愛好者加入：學習交流群 667279387