A basic problem in optics concerns lenses that map a source into an image. If we imagine the source image as a collection of point sources of light, a "perfect" lens would map each unique point in the source plane onto a unique point in the image plane. Ignoring diffraction effects, this is precisely what a pinhole camera does.
But real lenses are not perfect. Gravitational lenses, for instance, are pretty "blurry". Light from a point source is spread out over a broad area in the image plane. This behavior is characterized by the so-called point-spread function, or PSF, of the lens.
Quite often, the PSF is not known. But in an important subset of cases, the PSF is known exactly. It is, then, possible in principle to reconstruct the original source from a blurry image, but not without cost: The reconstruction can substantially increase noise. It is important to be able to estimate this "penalty".
Let's start with the basics. Let us characterize the source plane by $I({\bf x})$ representing the intensity of light at the location ${\bf x}$. Similarly, let us characterize the image plane by $I'({\bf x}')$. Let the PSF be ${\rm PSF}({\bf x},{\bf x}')$. That is, the PSF acts as a transfer function, characterizing how much light from a point source at ${\bf x}$ is deposited at location ${\bf x}'$ in the image plane. (Because ultimately, the PSF characterizes behavior associated with point sources, it is also sometimes called the impulse response function.)
The relationship between source and image is, then, given by
$$I'({\bf x}')=\iint d^2{\rm x}~I({\bf x}){\rm PSF}({\bf x},{\bf x}').\tag{1}$$
This is obviously the definition of an integral operator, the convolution operator:
$$I'({\bf x}')={\cal C}({\bf x},{\bf x}')\star I({\bf x}).\tag{2}$$
If we knew how to invert this operator, we'd know how to reconstruct $I({\bf x})$ from $I'({\bf x}')$:
$$I({\bf x})=\overline{{\cal C}({\bf x}',{\bf x})}\star I'({\bf x}').\tag{3}$$
To do this, we recognize that ${\cal C}$ is a linear operator. Discretizing the problem and representing the source and the image at $1\le i\le N$ and $1\le j\le M$ discrete locations, respectively, we can write the convolution as a sum:
$$I'({\bf x}'_j)=\sum\limits_{i=1}^N{\cal C}_{ij}I({\bf x}_i),\tag{4}$$
where
$${\cal C}_{ij}={\rm PSF}({\bf x}_i,{\bf x}_j').\tag{5}$$
This representation remains valid so long as the $N$ discrete locations fully and uniformly cover the entire source image; that is to say, they represent area elements of equal size, the area elements are sufficiently small such we do not need to worry about variations of intensity within an element, and no luminous part of the source remains uncovered. In short, we replaced the double integral over the entire source plane with a simple summation over the luminous parts, using a standard rectangle approximation of the integral.
If $N=M$, then ${\cal C}_{ij}$ is a square matrix. If that square matrix is degenerate, information is irretrievably lost and deconvolution is not possible. If the square matrix is invertible, then deconvolution is possible, and it is accomplished simply as
$$I({\bf x}_i)=\sum\limits_{j=1}^N\overline{{\cal C}_{ji}}I'({\bf x}_j').\tag{6}$$
In real-world problems, however, the image is not noise-free: $I'({\bf x}_j')\to I'({\bf x}_j')+\eta({\bf x_j'}),$ where $\eta$ represents random measurement noise. Usually, it is assumed that the noise is approximately uniform and follows a normal distribution. Furthermore, the noise is assumed to be uncorrelated: there is no correlation between $\eta({\bf x}_j')$ and $\eta({\bf x}_k')$ when $j\ne k$.
If these conditions are true, then as we apply the deconvolution operator, Gaussian noise must be added using root-mean-square summation:
$$I''({\bf x}_i)=\sum\limits_{j=1}^N\overline{{\cal C}_{ji}}I'({\bf x}_j')+\sqrt{\sum\limits_{j=1}^N\overline{{\cal C}_{ji}^2}\eta({\bf x}_j')^2}=\sum\limits_{j=1}^N\overline{{\cal C}_{ji}}I'({\bf x}_j')+\sqrt{\sum\limits_{j=1}^N\overline{{\cal C}_{ji}^2}}~~\eta.\tag{7}$$
What we would like to know is the change in the signal-to-noise ratio, or SNR. The SNR in this case is defined as the signal average divided by the standard deviation, which is to say, the RMS average of noise. Before deconvolution, we have
$${\rm SNR}'=\frac{\dfrac{1}{N}\sum\limits_{j=1}^NI'({\bf x}_j')}{\sqrt{\dfrac{1}{N}\sum\limits_{j=1}^N\eta({\bf x}_j')^2}}=\frac{\langle I'\rangle}{\eta}.\tag{8}$$
After deconvolution, we have
$${\rm SNR}''=\frac{\dfrac{1}{N}\sum\limits_{i=1}^NI({\bf x}_i)}{\sqrt{\dfrac{1}{N}\sum\limits_{i=1}^N\sum\limits_{j=1}^N\overline{{\cal C}_{ij}^2}\eta({\bf x}_j')^2}}=\frac{\langle I\rangle}{\sqrt{\dfrac{1}{N}\sum\limits_{i=1}^N\sum\limits_{j=1}^N\overline{{\cal C}_{ij}^2}}~~\eta}.\tag{9}$$
We define the deconvolution penalty as
$$\frac{{\rm SNR}''}{{\rm SNR}'}.\tag{10}$$
In many cases of interest, the PSF can be written in the form
$${\rm PSF}({\bf x}_i,{\bf x}_j')={\rm PSF}(|{\bf x}_i+\alpha{\bf x}'_j|)\tag{11}$$
with $\alpha$ constant, and the discretization defined such that ${\bf x}_i+\alpha{\bf x}_i'=0.$ That is to say, the PSF in this case depends only on the distance in the image plane between the projected location of the source and the location where the signal is measured. The limiting case of a pinhole camera is, in fact, given by ${\rm PSF}({\bf x}_i,{\bf x}_j')=\delta^2(|{\bf x}_i+\alpha{\bf x}'_j|)$ where $\delta$ is the Dirac delta function. This notation makes it apparent that $\alpha$ characterizes the geometric projection from the source plane to the image plane.
Moreover, in many cases of interest, the convolution matrix is dominated by its diagonal. In these cases, the convolution matrix can be conservatively approximated by the form
$${\cal C}_{ij}\sim\tilde{\cal C}_{ij}=\mu{\mathbb I}_{ij}+\nu{\mathbb U}_{ij},\tag{12}$$
where ${\mathbb I}$ is the identity (unit) matrix, ${\mathbb U}$ is the "everywhere one" matrix, $\mu\gg \nu$ and
$$\begin{align*}
\mu&={\cal C}_{ii}-\nu={\rm PSF}(0)-\nu,\tag{13}\\
\nu&=\frac{1}{N^2}\sum\limits_{i=1}^N\sum\limits_{j=1}^N{\cal C}_{ij}\sim\frac{1}{AA'}\iint_A d^2{\bf x}\iint_{A'} d^2{\bf x'}~{\rm PSF}({\bf x},{\bf x}'),\tag{14}\end{align*}$$
where $A$ and $A'$ represent finite integration areas in the source and image planes, respectively.
The inverse of $\tilde{\cal C}_{ij}$ can be computed easily:
$$\overline{\tilde{\cal C}_{ij}}=\frac{1}{\mu}{\mathbb I}_{ij}-\frac{\nu}{\mu(\mu+N\nu)}{\mathbb U}_{ij}.\tag{15}$$
Moreover,
$$\langle I'\rangle=\frac{1}{N}\sum_{j=1}^N\sum\limits_{i=1}^N\tilde{\cal C}_{ij}I({\bf x}_i)=N\nu\langle I\rangle,\tag{16}$$
and
$$\frac{1}{N}\sum\limits_{i=1}^N\sum\limits_{j=1}^N\overline{\tilde{\cal C}_{ij}^2}=\left[\frac{\mu+N\nu-\nu}{\mu(\mu+N\nu)}\right]^2+(N-1)\left[\frac{\nu}{\mu(\mu+N\nu)}\right]^2\sim\frac{1}{\mu^2}.\tag{17}$$
Therefore,
$${\rm SNR}''\sim\frac{\mu\langle I\rangle}{\eta}\tag{18}$$
and the deconvolution penalty is calculated as
$$\frac{{\rm SNR}''}{{\rm SNR}'}\sim\frac{\mu\langle I\rangle}{\eta}\frac{\eta}{N\nu\langle I\rangle}=\frac{\mu}{N\nu}.\tag{19}$$
In particular, if, for instance, $\mu\sim \sqrt{N}\nu$ (which is the case that appears to characterize gravitational lenses), we have
$$\frac{{\rm SNR}''}{{\rm SNR}'}\sim\frac{1}{\sqrt{N}}.\tag{20}$$
