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

Question

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?

how about instead of using the full Covar matrix 512^2, use the diagonal form, which is a vector of 512. — GameOfThrows
I didn't understand, can you elaborate? And also, how the noise will appear in image? — Ashutosh Gupta
your white noise vector is a matrix of 512^2 and I assume you are doing additive noise. Based on matlab specification "If A is a matrix whose columns represent random variables and whose rows represent observations, C is the covariance matrix with the corresponding column variances along the diagonal." Then cov(A) will simply return the covariance in the diagonal form. As for if the assumption stands, it is debatable. — GameOfThrows
Then how do I find entry C(k,l) of covariance between w_k and w_l ? My C matrix will be 512^2 X 512^2. — Ashutosh Gupta

kmac kmac · Accepted Answer · 2015-08-25T06:22:36

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;

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

1 Answers