# HG changeset patch # User Ivan Damnjanovic lnx # Date 1311864572 -3600 # Node ID 485747bf39e0de934f92586589897ac5bac3fcae # Parent fec205ec6ef699ecd2acfa2ef3ab92cb73ed13cc Two step dictonary learning - Integration of the code for dictionary update and dictionary decorrelation from Boris Mailhe diff -r fec205ec6ef6 -r 485747bf39e0 DL/two-step DL/SMALL_two_step_DL.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/two-step DL/SMALL_two_step_DL.m Thu Jul 28 15:49:32 2011 +0100 @@ -0,0 +1,127 @@ +function DL=SMALL_two_step_DL(Problem, DL) + +% determine which solver is used for sparse representation % + +solver = DL.param.solver; + +% determine which type of udate to use ('KSVD', 'MOD', 'ols' or 'mailhe') % + +typeUpdate = DL.name; + +sig = Problem.b; + +% determine dictionary size % + +if (isfield(DL.param,'initdict')) + if (any(size(DL.param.initdict)==1) && all(iswhole(DL.param.initdict(:)))) + dictsize = length(DL.param.initdict); + else + dictsize = size(DL.param.initdict,2); + end +end +if (isfield(DL.param,'dictsize')) % this superceedes the size determined by initdict + dictsize = DL.param.dictsize; +end + +if (size(sig,2) < dictsize) + error('Number of training signals is smaller than number of atoms to train'); +end + + +% initialize the dictionary % + +if (isfield(DL.param,'initdict')) + if (any(size(DL.param.initdict)==1) && all(iswhole(DL.param.initdict(:)))) + dico = sig(:,DL.param.initdict(1:dictsize)); + else + if (size(DL.param.initdict,1)~=size(sig,1) || size(DL.param.initdict,2) 1e-6); % ensure no zero data elements are chosen + perm = randperm(length(data_ids)); + dico = sig(:,data_ids(perm(1:dictsize))); +end + +% flow: 'sequential' or 'parallel'. If sequential, the residual is updated +% after each atom update. If parallel, the residual is only updated once +% the whole dictionary has been computed. Sequential works better, there +% may be no need to implement parallel. Not used with MOD. + +if isfield(DL.param,'flow') + flow = DL.param.flow; +else + flow = 'sequential'; +end + +% learningRate. If the type is 'ols', it is the descent step of +% the gradient (typical choice: 0.1). If the type is 'mailhe', the +% descent step is the optimal step*rho (typical choice: 1, although 2 +% or 3 seems to work better). Not used for MOD and KSVD. + +if isfield(DL.param,'learningRate') + learningRate = DL.param.learningRate; +else + learningRate = 0.1; +end + +% number of iterations (default is 40) % + +if isfield(DL.param,'iternum') + iternum = DL.param.iternum; +else + iternum = 40; +end +% determine if we should do decorrelation in every iteration % + +if isfield(DL.param,'coherence') + decorrelate = 1; + mu = DL.param.coherence; +else + decorrelate = 0; +end + +% show dictonary every specified number of iterations + +if (isfield(DL.param,'show_dict')) + show_dictionary=1; + show_iter=DL.param.show_dict; +else + show_dictionary=0; + show_iter=0; +end + +% This is a small patch that needs to be resolved in dictionary learning we +% want sparse representation of training set, and in Problem.b1 in this +% version of software we store the signal that needs to be represented +% (for example the whole image) + +tmpTraining = Problem.b1; +Problem.b1 = sig; +Problem = rmfield(Problem, 'reconstruct'); +solver.profile = 0; + +% main loop % + +for i = 1:iternum + solver = SMALL_solve(Problem, solver); + [dico, solver.solution] = dico_update(dico, sig, solver.solution, ... + typeUpdate, flow, learningRate); + if (decorrelate) + dico = dico_decorr(dico, mu, solver.solution); + end + Problem.A = dico; + if ((show_dictionary)&&(mod(i,show_iter)==0)) + dictimg = SMALL_showdict(dico,[8 8],... + round(sqrt(size(dico,2))),round(sqrt(size(dico,2))),'lines','highcontrast'); + figure(2); imagesc(dictimg);colormap(gray);axis off; axis image; + pause(0.02); + end +end + +Problem.b1 = tmpTraining; +DL.D = dico; + +end \ No newline at end of file diff -r fec205ec6ef6 -r 485747bf39e0 DL/two-step DL/dico_color.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/two-step DL/dico_color.m Thu Jul 28 15:49:32 2011 +0100 @@ -0,0 +1,46 @@ +function colors = dico_color(dico, mu) + % DICO_COLOR cluster a dictionary in pairs of high correlation atoms. + % Called by dico_decorr. + % + % Parameters: + % -dico: the dictionary + % -mu: the correlation threshold + % + % Result: + % -colors: a vector of indices. Two atoms with the same color have a + % correlation greater than mu + + numAtoms = length(dico); + colors = zeros(numAtoms, 1); + + % compute the correlations + G = abs(dico'*dico); + G = G-eye(size(G)); + + % iterate on the correlations higher than mu + c = 1; + maxCorr = max(max(G)); + while maxCorr > mu + % find the highest correlated pair + x = find(max(G)==maxCorr, 1); + y = find(G(x,:)==maxCorr, 1); + + % color them + colors(x) = c; + colors(y) = c; + c = c+1; + + % make sure these atoms never get selected again + G(x,:) = 0; + G(:,x) = 0; + G(y,:) = 0; + G(:,y) = 0; + + % find the next correlation + maxCorr = max(max(G)); + end + + % complete the coloring with singletons + index = find(colors==0); + colors(index) = c:c+length(index)-1; +end diff -r fec205ec6ef6 -r 485747bf39e0 DL/two-step DL/dico_decorr.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/two-step DL/dico_decorr.m Thu Jul 28 15:49:32 2011 +0100 @@ -0,0 +1,48 @@ +function dico = dico_decorr(dico, mu, amp) + %DICO_DECORR decorrelate a dictionary + % Parameters: + % dico: the dictionary + % mu: the coherence threshold + % amp: the amplitude coefficients, only used to decide which atom to + % project + % + % Result: + % dico: a dictionary close to the input one with coherence mu. + + % compute atom weights + if nargin > 2 + rank = sum(amp.*amp, 2); + else + rank = randperm(length(dico)); + end + + % several decorrelation iterations might be needed to reach global + % coherence mu. niter can be adjusted to needs. + niter = 1; + while niter < 5 && ... + max(max(abs(dico'*dico -eye(length(dico))))) > mu + 10^-6 + % find pairs of high correlation atoms + colors = dico_color(dico, mu); + + % iterate on all pairs + nbColors = max(colors); + for c = 1:nbColors + index = find(colors==c); + if numel(index) == 2 + % decide which atom to change (the one with lowest weight) + if rank(index(1)) < rank(index(2)) + index = fliplr(index); + end + + % update the atom + corr = dico(:,index(1))'*dico(:,index(2)); + alpha = sqrt((1-mu*mu)/(1-corr*corr)); + beta = corr*alpha-mu*sign(corr); + dico(:,index(2)) = alpha*dico(:,index(2))... + -beta*dico(:,index(1)); + end + end + niter = niter+1; + end +end + diff -r fec205ec6ef6 -r 485747bf39e0 DL/two-step DL/dico_update.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/two-step DL/dico_update.m Thu Jul 28 15:49:32 2011 +0100 @@ -0,0 +1,107 @@ +function [dico, amp] = dico_update(dico, sig, amp, type, flow, rho) + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % [dico, amp] = dico_update(dico, sig, amp, type, flow, rho) + % + % perform one iteration of dictionary update for dictionary learning + % + % parameters: + % - dico: the initial dictionary with atoms as columns + % - sig: the training data + % - amp: the amplitude coefficients as a sparse matrix + % - type: the algorithm can be one of the following + % - ols: fixed step gradient descent + % - mailhe: optimal step gradient descent (can be implemented as a + % default for ols?) + % - MOD: pseudo-inverse of the coefficients + % - KSVD: already implemented by Elad + % - flow: 'sequential' or 'parallel'. If sequential, the residual is + % updated after each atom update. If parallel, the residual is only + % updated once the whole dictionary has been computed. Sequential works + % better, there may be no need to implement parallel. Not used with + % MOD. + % - rho: learning rate. If the type is 'ols', it is the descent step of + % the gradient (typical choice: 0.1). If the type is 'mailhe', the + % descent step is the optimal step*rho (typical choice: 1, although 2 + % or 3 seems to work better). Not used for MOD and KSVD. + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + if ~ exist( 'rho', 'var' ) || isempty(rho) + rho = 0.1; + end + + if ~ exist( 'flow', 'var' ) || isempty(flow) + flow = sequential; + end + + res = sig - dico*amp; + nb_pattern = size(dico, 2); + + switch type + case 'rand' + x = rand(); + if x < 1/3 + type = 'MOD'; + elseif type < 2/3 + type = 'mailhe'; + else + type = 'KSVD'; + end + end + + switch type + case 'MOD' + G = amp*amp'; + dico2 = sig*amp'*G^-1; + for p = 1:nb_pattern + n = norm(dico2(:,p)); + % renormalize + if n > 0 + dico(:,p) = dico2(:,p)/n; + amp(p,:) = amp(p,:)*n; + end + end + case 'ols' + for p = 1:nb_pattern + grad = res*amp(p,:)'; + if norm(grad) > 0 + pat = dico(:,p) + rho*grad; + pat = pat/norm(pat); + if nargin >5 && strcmp(flow, 'sequential') + res = res + (dico(:,p)-pat)*amp(p,:); %#ok<*NASGU> + end + dico(:,p) = pat; + end + end + case 'mailhe' + for p = 1:nb_pattern + grad = res*amp(p,:)'; + if norm(grad) > 0 + pat = (amp(p,:)*amp(p,:)')*dico(:,p) + rho*grad; + pat = pat/norm(pat); + if nargin >5 && strcmp(flow, 'sequential') + res = res + (dico(:,p)-pat)*amp(p,:); + end + dico(:,p) = pat; + end + end + case 'KSVD' + for p = 1:nb_pattern + index = find(amp(p,:)~=0); + if ~isempty(index) + patch = res(:,index)+dico(:,p)*amp(p,index); + [U,S,V] = svd(patch); + if U(:,1)'*dico(:,p) > 0 + dico(:,p) = U(:,1); + else + dico(:,p) = -U(:,1); + end + dico(:,p) = dico(:,p)/norm(dico(:,p)); + amp(p,index) = dico(:,p)'*patch; + if nargin >5 && strcmp(flow, 'sequential') + res(:,index) = patch-dico(:,p)*amp(p,index); + end + end + end + end +end + diff -r fec205ec6ef6 -r 485747bf39e0 examples/Image Denoising/SMALL_ImgDenoise_DL_test_KSVDvsRLSDLAvsTwoStepMOD.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/examples/Image Denoising/SMALL_ImgDenoise_DL_test_KSVDvsRLSDLAvsTwoStepMOD.m Thu Jul 28 15:49:32 2011 +0100 @@ -0,0 +1,288 @@ +%% Dictionary Learning for Image Denoising - KSVD vs Recursive Least Squares +% +% This file contains an example of how SMALLbox can be used to test different +% dictionary learning techniques in Image Denoising problem. +% It calls generateImageDenoiseProblem that will let you to choose image, +% add noise and use noisy image to generate training set for dictionary +% learning. +% Two dictionary learning techniques were compared: +% - KSVD - M. Elad, R. Rubinstein, and M. Zibulevsky, "Efficient +% Implementation of the K-SVD Algorithm using Batch Orthogonal +% Matching Pursuit", Technical Report - CS, Technion, April 2008. +% - RLS-DLA - Skretting, K.; Engan, K.; , "Recursive Least Squares +% Dictionary Learning Algorithm," Signal Processing, IEEE Transactions on, +% vol.58, no.4, pp.2121-2130, April 2010 +% + + +% Centre for Digital Music, Queen Mary, University of London. +% This file copyright 2011 Ivan Damnjanovic. +% +% This program is free software; you can redistribute it and/or +% modify it under the terms of the GNU General Public License as +% published by the Free Software Foundation; either version 2 of the +% License, or (at your option) any later version. See the file +% COPYING included with this distribution for more information. +% +%% + + + +% If you want to load the image outside of generateImageDenoiseProblem +% function uncomment following lines. This can be useful if you want to +% denoise more then one image for example. +% Here we are loading test_image.mat that contains structure with 5 images : lena, +% barbara,boat, house and peppers. +clear; +TMPpath=pwd; +FS=filesep; +[pathstr1, name, ext, versn] = fileparts(which('SMALLboxSetup.m')); +cd([pathstr1,FS,'data',FS,'images']); +load('test_image.mat'); +cd(TMPpath); + +% Deffining the noise levels that we want to test + +noise_level=[10 20 25 50 100]; + +% Here we loop through different noise levels and images + +for noise_ind=2:2 +for im_num=1:1 + +% Defining Image Denoising Problem as Dictionary Learning +% Problem. As an input we set the number of training patches. + +SMALL.Problem = generateImageDenoiseProblem(test_image(im_num).i, 40000, '',256, noise_level(noise_ind)); +SMALL.Problem.name=int2str(im_num); + +Edata=sqrt(prod(SMALL.Problem.blocksize)) * SMALL.Problem.sigma * SMALL.Problem.gain; +maxatoms = floor(prod(SMALL.Problem.blocksize)/2); + +% results structure is to store all results + +results(noise_ind,im_num).noisy_psnr=SMALL.Problem.noisy_psnr; + +%% +% Use KSVD Dictionary Learning Algorithm to Learn overcomplete dictionary + +% Initialising Dictionary structure +% Setting Dictionary structure fields (toolbox, name, param, D and time) +% to zero values + +SMALL.DL(1)=SMALL_init_DL(); + +% Defining the parameters needed for dictionary learning + +SMALL.DL(1).toolbox = 'KSVD'; +SMALL.DL(1).name = 'ksvd'; + +% Defining the parameters for KSVD +% In this example we are learning 256 atoms in 20 iterations, so that +% every patch in the training set can be represented with target error in +% L2-norm (Edata) +% Type help ksvd in MATLAB prompt for more options. + + +SMALL.DL(1).param=struct(... + 'Edata', Edata,... + 'initdict', SMALL.Problem.initdict,... + 'dictsize', SMALL.Problem.p,... + 'exact', 1, ... + 'iternum', 20,... + 'memusage', 'high'); + +% Learn the dictionary + +SMALL.DL(1) = SMALL_learn(SMALL.Problem, SMALL.DL(1)); + +% Set SMALL.Problem.A dictionary +% (backward compatiblity with SPARCO: solver structure communicate +% only with Problem structure, ie no direct communication between DL and +% solver structures) + +SMALL.Problem.A = SMALL.DL(1).D; +SMALL.Problem.reconstruct = @(x) ImgDenoise_reconstruct(x, SMALL.Problem); + +%% +% Initialising solver structure +% Setting solver structure fields (toolbox, name, param, solution, +% reconstructed and time) to zero values + +SMALL.solver(1)=SMALL_init_solver; + +% Defining the parameters needed for image denoising + +SMALL.solver(1).toolbox='ompbox'; +SMALL.solver(1).name='omp2'; +SMALL.solver(1).param=struct(... + 'epsilon',Edata,... + 'maxatoms', maxatoms); + +% Denoising the image - find the sparse solution in the learned +% dictionary for all patches in the image and the end it uses +% reconstruction function to reconstruct the patches and put them into a +% denoised image + +SMALL.solver(1)=SMALL_solve(SMALL.Problem, SMALL.solver(1)); + +% Show PSNR after reconstruction + +SMALL.solver(1).reconstructed.psnr + +%% +% For comparison purposes we will denoise image with overcomplete DCT +% here +% Set SMALL.Problem.A dictionary to be oDCT (i.e. Problem.initdict - +% since initial dictionaruy is already set to be oDCT when generating the +% denoising problem + + +% Initialising solver structure +% Setting solver structure fields (toolbox, name, param, solution, +% reconstructed and time) to zero values + +SMALL.solver(2)=SMALL_init_solver; + +% Defining the parameters needed for image denoising + +SMALL.solver(2).toolbox='ompbox'; +SMALL.solver(2).name='omp2'; +SMALL.solver(2).param=struct(... + 'epsilon',Edata,... + 'maxatoms', maxatoms); + +% Initialising Dictionary structure +% Setting Dictionary structure fields (toolbox, name, param, D and time) +% to zero values + +SMALL.DL(2)=SMALL_init_DL('TwoStepDL', 'MOD', '', 1); + + +% Defining the parameters for MOD +% In this example we are learning 256 atoms in 20 iterations, so that +% every patch in the training set can be represented with target error in +% L2-norm (EData) +% Type help ksvd in MATLAB prompt for more options. + + +SMALL.DL(2).param=struct(... + 'solver', SMALL.solver(2),... + 'initdict', SMALL.Problem.initdict,... + 'dictsize', SMALL.Problem.p,... + 'iternum', 40,... + 'mu', 0.7,... + 'show_dict', 1); + +% Learn the dictionary + +SMALL.DL(2) = SMALL_learn(SMALL.Problem, SMALL.DL(2)); + +% Set SMALL.Problem.A dictionary +% (backward compatiblity with SPARCO: solver structure communicate +% only with Problem structure, ie no direct communication between DL and +% solver structures) + +SMALL.Problem.A = SMALL.DL(2).D; +SMALL.Problem.reconstruct = @(x) ImgDenoise_reconstruct(x, SMALL.Problem); + +% Denoising the image - find the sparse solution in the learned +% dictionary for all patches in the image and the end it uses +% reconstruction function to reconstruct the patches and put them into a +% denoised image + +SMALL.solver(2)=SMALL_solve(SMALL.Problem, SMALL.solver(2)); + +%% +% In the b1 field all patches from the image are stored. For RLS-DLA we +% will first exclude all the patches that have l2 norm smaller then +% threshold and then take min(40000, number_of_remaining_patches) in +% ascending order as our training set (SMALL.Problem.b) + +X=SMALL.Problem.b1; +X_norm=sqrt(sum(X.^2, 1)); +[X_norm_sort, p]=sort(X_norm); +p1=p(X_norm_sort>Edata); +if size(p1,2)>40000 + p2 = randperm(size(p1,2)); + p2=sort(p2(1:40000)); + size(p2,2) + SMALL.Problem.b=X(:,p1(p2)); +else + size(p1,2) + SMALL.Problem.b=X(:,p1); + +end + +% Forgetting factor for RLS-DLA algorithm, in this case we are using +% fixed value + +lambda=0.9998 + +% Use Recursive Least Squares +% to Learn overcomplete dictionary + +% Initialising Dictionary structure +% Setting Dictionary structure fields (toolbox, name, param, D and time) +% to zero values + +SMALL.DL(3)=SMALL_init_DL(); + +% Defining fields needed for dictionary learning + +SMALL.DL(3).toolbox = 'SMALL'; +SMALL.DL(3).name = 'SMALL_rlsdla'; +SMALL.DL(3).param=struct(... + 'Edata', Edata,... + 'initdict', SMALL.Problem.initdict,... + 'dictsize', SMALL.Problem.p,... + 'forgettingMode', 'FIX',... + 'forgettingFactor', lambda,... + 'show_dict', 1000); + + +SMALL.DL(3) = SMALL_learn(SMALL.Problem, SMALL.DL(3)); + +% Initialising solver structure +% Setting solver structure fields (toolbox, name, param, solution, +% reconstructed and time) to zero values + +SMALL.Problem.A = SMALL.DL(3).D; +SMALL.Problem.reconstruct = @(x) ImgDenoise_reconstruct(x, SMALL.Problem); + +SMALL.solver(3)=SMALL_init_solver; + +% Defining the parameters needed for image denoising + +SMALL.solver(3).toolbox='ompbox'; +SMALL.solver(3).name='omp2'; +SMALL.solver(3).param=struct(... + 'epsilon',Edata,... + 'maxatoms', maxatoms); + + +SMALL.solver(3)=SMALL_solve(SMALL.Problem, SMALL.solver(3)); + +SMALL.solver(3).reconstructed.psnr + + +% show results % + +SMALL_ImgDeNoiseResult(SMALL); + +results(noise_ind,im_num).psnr.ksvd=SMALL.solver(1).reconstructed.psnr; +results(noise_ind,im_num).psnr.odct=SMALL.solver(2).reconstructed.psnr; +results(noise_ind,im_num).psnr.rlsdla=SMALL.solver(3).reconstructed.psnr; +results(noise_ind,im_num).vmrse.ksvd=SMALL.solver(1).reconstructed.vmrse; +results(noise_ind,im_num).vmrse.odct=SMALL.solver(2).reconstructed.vmrse; +results(noise_ind,im_num).vmrse.rlsdla=SMALL.solver(3).reconstructed.vmrse; +results(noise_ind,im_num).ssim.ksvd=SMALL.solver(1).reconstructed.ssim; +results(noise_ind,im_num).ssim.odct=SMALL.solver(2).reconstructed.ssim; +results(noise_ind,im_num).ssim.rlsdla=SMALL.solver(3).reconstructed.ssim; + +results(noise_ind,im_num).time.ksvd=SMALL.solver(1).time+SMALL.DL(1).time; +results(noise_ind,im_num).time.rlsdla.time=SMALL.solver(3).time+SMALL.DL(3).time; +clear SMALL; +end +end +% save results.mat results diff -r fec205ec6ef6 -r 485747bf39e0 util/SMALL_init_DL.m --- a/util/SMALL_init_DL.m Tue Jul 26 16:01:17 2011 +0100 +++ b/util/SMALL_init_DL.m Thu Jul 28 15:49:32 2011 +0100 @@ -1,4 +1,4 @@ -function DL = SMALL_init_DL(varargin) +function DL = SMALL_init_DL(toolbox, name, param, profile) %% Function initialise SMALL structure for Dictionary Learning. % Optional input variables: % toolbox - name of Dictionary Learning toolbox you want to use @@ -17,9 +17,27 @@ % %% -DL.toolbox=[]; -DL.name=[]; -DL.param=[]; +if ~ exist( 'toolbox', 'var' ) || isempty(toolbox) + DL.toolbox = []; +else + DL.toolbox = toolbox; +end +if ~ exist( 'name', 'var' ) || isempty(name) + DL.name = []; +else + DL.name = name; +end +if ~ exist( 'param', 'var' ) || isempty(param) + DL.param = []; +else + DL.param = param; +end +if ~ exist( 'profile', 'var' ) || isempty(profile) + DL.profile = 1; +else + DL.profile = profile; +end + DL.D=[]; DL.time=[]; end \ No newline at end of file diff -r fec205ec6ef6 -r 485747bf39e0 util/SMALL_learn.m --- a/util/SMALL_learn.m Tue Jul 26 16:01:17 2011 +0100 +++ b/util/SMALL_learn.m Thu Jul 28 15:49:32 2011 +0100 @@ -18,8 +18,9 @@ % License, or (at your option) any later version. See the file % COPYING included with this distribution for more information. %% - +if (DL.profile) fprintf('\nStarting Dictionary Learning %s... \n', DL.name); +end start=cputime; tStart=tic; if strcmpi(DL.toolbox,'KSVD') @@ -58,6 +59,17 @@ D(:,i)=D(:,i)/norm(D(:,i)); end + elseif strcmpi(DL.toolbox,'TwoStepDL') + + DL=SMALL_two_step_DL(Problem, DL); + + % we need to make sure that columns are normalised to + % unit lenght. + + for i = 1: size(DL.D,2) + DL.D(:,i)=DL.D(:,i)/norm(DL.D(:,i)); + end + D = DL.D; % To introduce new dictionary learning technique put the files in % your Matlab path. Next, unique name for your toolbox needs % to be defined and also prefferd API for toolbox functions @@ -83,9 +95,11 @@ %% % Dictionary Learning time tElapsed=toc(tStart); - DL.time = cputime - start; +DL.time = cputime - start; +if (DL.profile) fprintf('\n%s finished task in %2f seconds (cpu time). \n', DL.name, DL.time); fprintf('\n%s finished task in %2f seconds (tic-toc time). \n', DL.name, tElapsed); +end DL.time=tElapsed; % If dictionary is given as a sparse matrix change it to full