Image de-blurring: in frequency domain or in spatial domain ? 

Very often image deblurring, i.e. reconstruction ("sharpening")  of blurred images is performed in discrete frequency domain instead of spatial domain. Despite having many positive aspects,  deblurring in frequency domain has a negative feature of critical importance: the blurring matrix in the frequency domain is very ill-conditioned.

But what exactly “ill-conditioned” means, and what effect does it have on image restoration?

Consider in more detail the matrix condition ratio. It is defined as the ratio between the matrix' maximum and minimum singular values:

cond(A) = s_max/s_min

where s_max and s_min are the maximum and the minimum diagonal values of the matrix S in the singular value decomposition of our matrix A:

A = U S V^T

The meaning of the condition ratio should be considered in more detail. Consider a system of linear algebraic equations (or SLAE for short):  Ax=b.

Suppose we know the elements of the matrix A exactly, but elements of the right-hand side vector b are known with some degree of uncertainty b+db. The uncertainty in the elements of the right-hand side vector b leads to an uncertainty in the solution x. Then we have

A(x+dx)= b+db

Then, following Forsythe and Moler, we have

which means that  “uncertainty” in the solution of SLAE is limited by cond(A) multiplied by the “uncertainty” in the right-hand side vector b.

Similarly,

which means that  “uncertainty” in the solution of SLAE is limited by cond(A) multiplied by the “uncertainty” in the elements of the matrix A.

In practice, the above expressions mean that if cond(A) = 10⁶ then any “uncertainties” (or numerical discrepancies, to be more precise) in the elements of the right-hand side vector b or of the matrix A are multiplied in the solution by a value of the magnitude of about one million. 

The larger the condition ratio of the matrix is, the less accurate is the solution, and vice versa.

For some filters in the frequency domain, a threshold value is introduced. All elements of the blurring matrix, which are less or equal to that threshold values, are set to zero. This procedure effectively reduces the rank of the blurring matrix and improves its condition number. However, because small singular values of the blurring matrix correspond to high frequencies, which represent the details of the image, the fine details if the image are irreparably lost.

As an example, we consider the condition ratio of a simple blurring matrix in both spatial and frequency domains.

Consider an LSI degradation system with a one-dimensional motion blur in horizontal direction spread over 9 pixels and normalized.

The following MATLAB code reads an original RGB image of a rose from a .jpeg file, converts it to grey scale, and introduces blurring in frequency domain:

original = im2double(rgb2gray((imread('input_tiny.jpg'))));

% 1-D motion blur

motion_kernel = ones(1, 9) / 9;

% frequency response

motion_freq = fft2(motion_kernel, 64, 64);

cond = cond(motion_freq);

As expected, the condition ratio of the blurring matrix is very large:  8.9*10⁺¹⁸, which MATLAB sets to infinity, thus indicating that the matrix is singular and, therefore, no reliable solution could be obtained.

The condition ratio of the blurring matrix in the spatial domain (constructed when rows of the blurred image are ordered lexicographically) equals 36.686, which means that the matrix is very well-conditioned, and the solution will have sufficient (in fact, 13 or 14) significant digits. 

The  blurred and restored images are shown below:

The general conclusion is the image de-blurring in spatial domain is the only reliable approach to the problem. Working in frequency domain should be avoided at all costs. 

It should be noted that it does not matter what mathematical method is used to solve the SLAE with ill-conditioned matrix, since all the errors (uncertainties) in the input data as well as the round-off errors would be magnified in the solution.

An interesting linear algebra problem from signal processing. 

This problem was presented on https://stackoverflow.com/questions/52373413/solve-underdetermined-system-of-equations-for-sparse-solution/52404876#52404876 
It appears that a discussion devoted to it might be interesting, in order to analyse possible approaches to solution. 

We approach the problem by transforming the equation into a form more "recognizable" from the classic linear algebra approach to the solution of systems of linear equations. 

The last expression, essentially, represents a set of m systems of linear equations with one and the same matrix B^T and different right-hand sides, which are the columns of C^T. In order to proceed further, we need to perform Gaussian elimination for the matrix B^T , simultaneously with all the right-hand sides, using some pivoting algorithm. As the matrix is singular, the resulting system of equations, for each particular column j from C^T and, correspondingly, for each particular column j from A^T, will have the following form: 

Since the matrix of the system is singular, k<n , where k=rank(B^T) . The vector of unknowns consists of two smaller vectors: 

 of the size k , containing the components from 1st to k-th of the column j of the A^T

                matrix;

  of the size n-k, containing the components from (k+1)th to n-th of the same                           column of  A^T. 

Correspondingly:

  of the size k, containing (permuted, due to the Gaussian elimination with column                    pivoting performed) components of the j-th column of C^T ;

  of the size n-k, containing the remaining components  of the j-th column of C^T . 

Therefore, upon completion of the Gaussian elimination, the last few equations in the system are trivial:

Therefore, two possibilities only exist: 


1) All of the components of the lower part of the right-hand side vector are zeros:

In that case, the system has infinite number of solutions, since any alpha-vector would satisfy the subsystem. Since our goal is to find a solution of the original system as "sparse" as possible, we can safely set all components of the alpha-vector to zeros, and solve the remaining equations by backward substitution, thus obtaining a unique solution. 

2) At least one component of the right-hand side (ie. of the chi-vector) of the above subsystem is different from zero. In that case, the subsystem has no solutions at all, and, therefore, the original system has no solutions as well. In such case, the only way to proceed is to try to find an approximation to the solution in terms of the least squares. 

We will perform the analysis from the very beginning, using a well-established method of Singular Value Decomposition, or SVD for short. Given a system of linear equations, say, Fx=g , matrix F can be represented as the product of multiplication of three matrices (we limit ourselves by case of square matrix F , i.e. in the context of the current problem): 

In the context of the problem which is being discussed here, the diagonal matrix does contain at least one zero element on the diagonal, since the original matrix is supposed to be singular. Here, we will consider a simple SLAE of the size 5x5 where 

- two rows of the matrix are linearly dependent, and

- the system is incompatible, so that no exact solution exists, and therefore we will try to find a least squares approximation to the solution.

The 5-th row of the matrix equals the product of multiplication of the 3rd row by the factor of -2. Now, if we select the right-hand side of the system in such a way that the 5th element of the vector g does not equal the 3rd element multiplied by -2 , the whole SLAE is incompatible, and we can only try to find a least squares approximation to the solution. 

Lets put the right-hand side vector g as follows: g = [ 1 ,  2,  3,  4,  5 ]^T 

The SLAE is, obviously, incompatible. 

Now, what the SVD algorithm delivers:

The elements of the diagonal matrix, named here as singular values, contain one zero element, as expected. The SVD-algorithm prescribes the following procedure to obtain a least-squares solution of the system:

1) Introduce vector d=U^Tg .

2) Calculate vector z as follows:

   

                                z_j = d_j / s_j  , if s_j is a non-zero value;

                                z_j can be chosen arbitrarily, if s_j = 0. 

3) The solution is obtained as the following product: x=Vz . 

And here comes the beauty of the SVD-method: regardless of the chosen value of z_j for which  s_j = 0, the program will deliver a solution vector from the R⁵ space such that the norm of the residuals vector will be the same, and it will be the least possible value for this system. (Obviously, if the system was compatible, the residuals would all equal zero). 

We provide just three cases to demonstrate that for our example. Namely, we set z₅ to 0., to 1., and to a randomly selected value, and see what happens. 

As we can see, each particular choice of z₅ leads to a different solution vector. Nevertheless, the residuals vector is one and the same in all cases, which means that all three solutions are equally valid in terms of least squares. 

Since there is an infinite number of least-squares solutions, each particular one corresponding to one particular value of chosen z₅, we can probably consider an opposite task: set one component of the solution vector to zero, and find another components such that the vector of the residuals will stay the same (obviously, some particular value of z will be calculated along the way). The more singular values s_i  equal zero, the more components of the solution vector can theoretically be made equal zero. We should recall that our original goal was to deliver as sparse matrix A as possible. We consider that as possible subject of further algorithm development. 

                       

Solution of least squares problems using normal equations. 

Least squares problems assume that an over-determined system of linear algebraic equations needs to be solved in a way that the solution vector somehow minimizes the residuals vector that arises when the actual right-hand side of the system is subtracted from the product of the system's matrix and the obtained solution. 

A system of linear algebraic equations is over-determined if the number of rows in the system exceeds the number of columns (a case when some of the equations are linearly dependent is omitted for simplicity). Solving such system in classical sense is impossible; however, a method is known that transforms such system into a system of equations with square, symmetric, positive definite matrix. The obtained system is called a system of normal equations. The procedure is as follows: the original system is multiplied from the left by the original matrix transpose. 

The resulting matrix of the system, since it is now symmetric and positive definite, allows application of the Conjugate Gradients Method (or CGM). That fact has a very big significance for large-scale problems with sparse matrices: the CGM allows solution of systems consisting of millions of equations, and the matrix can be stored completely in the computer's Random Access Memory (or RAM) as long as the matrix is sparse (ie. if it contains, say, a few dozens of non-zero elements per row). 

In order to illustrate how the CGM can be applied to over-determined matrices, and how the numerical quality of the solution can be evaluated, we consider a rather simple example.

In order to obtain a solution, the approach based on the use of normal equations  and the Conjugate Gradients Method has been used. The solver was written in C and run on a computer with Intel Celeron N3450 CPU. 

An issue of special importance is verification of the solution. Even if the CGM iterations have converged, we can have reasonable doubts whether the obtained solution is numerically solid. The way to perform such verification is to multiply the original matrix by the solution vector and check: how different the obtained right-hand side (or RHS) is from the original one. As it can be seen from the example above, the solution is trustworthy, since the b[i] vector components aren't different from the original ones.