CODEWITHSUNDEEP
matlabnoise
Ashutosh Gupta
Ashutosh Gupta

Reputation: 680

How to add general [with a given covariance matrix] noise in images in MATLAB?

I am trying to work with additive non-white model for denoising. The problem is that I need full covariance matrix of the noise.

I want to add non-white noise to my 512 X 512 (or larger) image. The white noise vector length will be 512^2. Hence the covariance matrix will be prohibitively large (512^2 X 512^2). I cannot supply covar matrix to mvrnd to get a vector since I cannot even store it. My question is:

How can I add non-white noise while getting an expression for the j,k entry of my covariance matrix as a formula?

Upvotes: 2

Views: 1383

Answers (1)

kmac
kmac

Reputation: 698

This problem interested me, so I dug in. If you can use a sparse inverse covariance matrix for your noise (note: the covariance matrix can still be dense), you can avoid ever storing the full covariance in memory. It just takes some tricks with sparse solvers and the Cholesky decomposition. See below:

Output with n = 64^2:

Timing for 4096x4096: sparse: 0.73495 ms, dense: 8.2234 ms, total: 3103.46 ms, error: 1.95399e-14 nnz in Cinv: 16268 in chol(Cinv): 603667

Name           Size                  Bytes  Class     Attributes

CLfull      4096x4096            134217728  double              
CinvL       4096x4096             11593432  double    sparse

Output with n = 512^2:

Timing for 262144x262144: sparse: 3.84335 ms, dense: NaN ms, total: 195.025 ms, error: NaN nnz in Cinv: 409438 in chol(Cinv): 347078

Name            Size                  Bytes  Class     Attributes

CinvL      262144x262144            7703256  double    sparse

Code:

clear;
rng(0);
tic_total = tic;

%% Generate a random sparse inverse covariance matrix

n = 64^2; % 4096
%n = 256^2; % 65536
%n = 512^2; % 262144

% We limit the density to 1/n + 65536/(n^2) for large matrices
% passing positive singular values to sprand, and ensuring positive
% semi-definite by multiplying by its transpose
%Cinv = sprand(n,n,1/n+min(1/n,65536/n^2),rand(1,n)+1); % <-- This takes too long for large n
Cinv = sprand(n,n,min(1/n,65536/n^2)) + speye(n); % <-- Instead, add diagonal
Cinv = Cinv*Cinv'; %'

% Set mean to 0 (or something) ...
mu = zeros(n,1); % randn(n,1);

%% Cholesky decomposition of inverse covariance

CinvL = chol(Cinv, 'lower');
tc = toc;

%% Use a sparse solver to un-whiten the normal sample
% Drawing values from a normal distribution:
% https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Drawing_values_from_the_distribution
% Cholesky decomposition for (inverse) covariance matrix
% http://scicomp.stackexchange.com/questions/3188/dealing-with-the-inverse-of-a-positive-definite-symmetric-covariance-matrix

% Normal white sample
w = randn(n,1);

%% Sparse version
% Unwhitening is normal s = mu + CL*w, where CL = chol(C,'lower')
% But we have CinvL = chol(Cinv, 'lower') = inv(CL)
repeat = 20;
tic;
for i=1:repeat % Loop for average timing / caching funniness
  s1 = mu + CinvL\w; % <-- Sparse solver, fast and small memory footprint
end
t1 = toc/repeat;

%% Dense version
if n <= 4096 % Any larger and we might run out of memory
  CLfull = CinvL\eye(n);
  tic;
  for i=1:repeat
    s2 = mu + CLfull*w; % <-- Dense solver, slow with large memory footprint
  end
  t2 = toc/repeat;
else
  t2 = NaN;
  s2 = NaN*s1;
end

% Display timing and sizes
fprintf('Timing for %ux%u: sparse: %g ms, dense: %g ms, total: %g ms, error: %g\n', ...
  n, n, 1000*t1, 1000*t2, 1000*toc(tic_total), max(s1-s2));
fprintf('nnz in Cinv: %u in chol(Cinv): %u\n', nnz(Cinv), nnz(CinvL));
whos CinvL CLfull;

Upvotes: 1

Related Questions