Artur Sowa

Image processing via simulated quantum dynamics

De Gruyter | Published online: January 28, 2017


We apply simulated quantum evolution to image processing, and examine its practicality in the context of image denoising. More specifically, in our approach image processing consists of three stages: First, a digitized gray-scale image is represented as a quantum variable – typically, a density matrix. Second, the quantum variable is evolved via the Markovian master equation in Lindblad form. Third, the quantum variable is back-converted into an image. Numerical experiments indicate remarkable denoising results are obtained in this way for a suitable choice of flow parameters. To our knowledge the proposed image processing technique is conceptually new.

MSC 2010: 60G35; 81S22; 93E11

1 Introduction

Image processing (IP) literature abounds with applications of mathematical techniques of most sophisticated and innovative variety, see e.g. [5]. It has traditionally relied on techniques from Harmonic Analysis, Partial Differential Equations, Differential Geometry, Statistical Physics, Computer Optimization, and, more recently, reached for inspiration to Statistical Physics (see e.g. [10]) or Artificial Intelligence (see e.g. [11]). At the same time, there appear to have been virtually no attempts at engaging in IP bona fide quantum mechanical models. We believe the reason may have been, first and foremost, the counterintuitive structure of quantum mechanics itself. However, there are many incentives nowadays to give such an idea a thorough consideration. At the time when quantum technologies draw nigh and, in particular, quantum information processing is becoming a reality – e.g. with the availability of D-Wave Two quantum annealer [3] – a quantum-vantage point at IP will soon become a matter of necessity. It is therefore timely to explore the possibility of engaging the quantum point of view toward a solution of the traditional and fundamental IP challenges. The work presented in this paper reports some observations of relevance. They grew via serendipity out of the basic quest to interpret the effects of quantum dynamics in a classical setting, such as endowed by a natural image. To this end, we have utilized preexisting algorithms first proposed in [17]. In order to make this interdisciplinary material reasonably self-contained, we include a brief outline of the theoretical background that is necessary to rigorously explain our method. We then describe and visualize the results of numerical experiments. Our goal is to share with the readership the essential concepts rather than technical details, and so we defer the quantitative analysis of image enhancement parameters, which more appropriately belongs in engineering literature, to future communications. It is our hope that the proposed approach to IP contributes in two ways: first, it yields a novel classical IP technique based on simulated quantum dynamics; second, it possibly opens a way for IP based on a hard-wired quantum dynamics. With regards to the latter goal, it remains to be seen if science will find a way to implement data such as a natural image via a state of a quantum system, or a quantum observable. It seems natural to look for opportunities to carry out such a program first in those situations where the image naturally emerges from a quantum environment as it happens in the magnetic resonance imaging (MRI) or in the magnetoencephalography (MEG) machines. Here we only outline some of the mathematical realities that such a task will entail.

2 Representing images via quantum states

In this section we propose several ways of representing a digital gray-scale image as a quantum variable. The list is not exhaustive and more methods may be proposed. However, the numerical implementation of either one of the schemas proposed here requires only the Fourier transforms and rather trivial symmetrization. Thus, critically, these schemas may be implemented via fast and memory efficient algorithms. In a generalized scenario, not considered here, the Fourier transform may be replaced by other fast transforms, as long as certain conditions of symmetry are satisfied – see (1) below – such transforms have been described in [15, 14].

2.1 Representing a digital gray-scale image as a density matrix

Let us fix the standard Hermitian form in N . Namely, let

𝐯 | 𝐰 = k v ¯ k w k

for any two vectors 𝐯 , 𝐰 N , where 𝐯 = k v k 𝐞 k and { 𝐞 k } is the standard basis, etc. Recall that ρ N × N is a density matrix if it satisfies the following three conditions:

  1. (1)

    ρ is Hermitian, i.e. 𝐯 | ρ 𝐰 = ρ 𝐯 | 𝐰 for all 𝐯 , 𝐰 N ; briefly, ρ = ρ .

  2. (2)

    ρ is positive semidefinite, i.e. 𝐯 | ρ 𝐯 0 for all 𝐯 N .

  3. (3)

    ρ has trace 1, i.e. Tr ρ := ρ 1 , 1 + ρ 2 , 2 + + ρ N , N = 1 .

Next, let I be an N × N matrix with real entries, e.g. the luminance map of a digitized image. Our goal is to encode I in a density matrix, say, ρ N × N . Let I ^ be the discrete Fourier transform of I. Since I is real, I ^ has a point-symmetry, i.e.

(1) I ^ = X + rot 180 X + D ,

where D is the diagonal of I ^ , X is the upper triangular matrix containing the above-diagonal entries of I ^ , and rot 180 X denotes rotation of X about its center via 180 . Next, define

(2) ρ ~ = X + X + α D ˇ and ρ = β ρ ~ , where β = 1 Tr ρ ~ .

Here, D ˇ denotes the (one-dimensional) discrete inverse Fourier transform of D. In particular, D ˇ is real. The constant α > 0 is chosen in such a way as to ensure that ρ is positive semi-definite. Indeed, for α large enough ρ ~ is positive semi-definite – see the discussion in the Appendix (Section A). Thus, ρ defined in this way is a density matrix.

Note also that the operation I ρ – as defined via (1) and (2) – is invertible, provided one keeps in memory the additional parameters α and β. We will adopt the following shorthand notation for the I ρ function and its inverse:

(3) [ ρ , α , β ] = image2rho ( I , α ) , I = rho2image ( ρ , α , β ) .

Remark 2.1

Note that the construction of ρ given above cannot be generalized to the case of continuous variable I : [ 0 , 2 π ] × 2 . Indeed, in the latter case D ˇ : [ 0 , 2 π ] is not suitable for a matrix diagonal. In the next subsection we demonstrate that an image may be encoded in a pair of quantum observables in a way that admits generalization to the continuous variable case.

2.2 Representing an image via a pair of quantum observables

A quantum-mechanical observable is a self-adjoint operator which, unlike the state variable, need not be positive definite nor of trace one. Here, we propose a method of representing an image via a pair of self-adjoint operators. We will discuss the construction in the continuous setting I : [ 0 , 2 π ] × 2 , the natural reduction to the discrete case being straightforward. For simplicity we assume that I is square integrable.


Let L 2 [ 0 , 2 π ] × 2 be the Hilbert space of square integrable complex functions. We fix an orthonormal basis { ϕ k , l } , a modification of the Fourier basis, defined as follows:[1]

ϕ k , l ( x , y ) e i k x e - i l y , k , l .

The minus sign in the second exponent is intentional – it ensures the following crucial property:

(4) ϕ k , l ( x , y ) = ϕ l , k ( y , x ) ¯ .

For a function g L 2 [ 0 , 2 π ] × 2 we denote its basis coefficients z k , l [ g ] , i.e.

z k , l [ g ] = g | ϕ k , l = 0 2 π 0 2 π g ( x , y ) ϕ k , l ( x , y ) ¯ 𝑑 x 𝑑 y ,

so that

g ( x , y ) = k , l z k , l [ g ] ϕ k , l ( x , y ) .

Next, let f L 2 [ 0 , 2 π ] × 2 be a real-valued function satisfying the constraint

(5) f ( x , y ) = f ( y , x ) .

It is straightforward to see that the coefficients z k , l = z k , l [ f ] form a self-adjoint matrix, i.e. z k , l = z ¯ l , k .

Luminance as a pair of observables.

With this understood, we define f 1 , f 2 : [ 0 , 2 π ] × 2 as follows:

f 1 ( x , y ) = I ( x , y ) if x y , and subsequently f 1 ( x , y ) = f 1 ( y , x ) for x > y ,


f 2 ( x , y ) = I ( x , y ) if x > y , and subsequently f 2 ( x , y ) = f 2 ( y , x ) for x < y .

For completeness we also set f 2 ( x , x ) = 0 for all x . Note that f 1 agrees with I below the diagonal, f 2 agrees with I above the diagonal, and both satisfy (5). It is also clear that I can be reconstructed from the pair f 1 , f 2 . Finally, we define the pair of observables ( F 1 , F 2 ) as the operators whose matrix representations in the basis { ϕ k , l } , are given by matrices

F 1 = ( z k , l [ f 1 ] ) k , l , F 2 = ( z k , l [ f 2 ] ) k , l .

Clearly, f 1 (resp. f 2 ) can be reconstructed from F 1 (resp. F 2 ). We will adopt the following shorthand notation for the I ( F 1 , F 2 ) function and its inverse:

(6) [ F 1 , F 2 ] = image2observables ( I ) , I = observables2image ( F 1 , F 2 ) .

3 Quantum dynamics as a means for evolving the luminance function

As demonstrated above, the luminance function may be encoded in a density matrix ρ or in a pair of quantum observables. Quantum mechanics provides precise rules for the time evolution of these types of variables. Our overall goal is to understand if quantum-mechanical evolution can effect changes that are desirable from the IP point of view. It is not easy to resolve the question completely. That is because a quantum system may evolve in a plethora of different ways. The brief synopsis below is meant to delineate what can be expected. A reader well familiar with the formalism of Quantum Theory need only skim over this material.

3.1 The rudimentary quantum dynamics: An unlikely IP tool?

According to the postulates of Quantum Mechanics every closed (meaning isolated, i.e. not interacting with anything else) quantum system is described via a Hilbert space, say, that encodes the totality of the system’s degrees of freedom. Furthermore, in a non-relativistic setting, which is the only one considered here, the time variable t is extrinsic. Every quantum system evolves in time continually changing its state. The state is mathematically expressed as a single (norm one) time-dependent vector ψ = ψ ( t ) (typically dressed as the so-called ket | ψ ). The time evolution of | ψ depends on the Hamiltonian, denoted H, which is either a bounded or an unbounded, but always self-adjoint, operator in . The evolution is described by the Schrödinger equation, which (in conveniently chosen physical units) assumes the form

i | ψ ˙ = H | ψ .

Here the dot above a variable signifies the time derivative. This equation is the single most fundamental object for the so-called circuit model quantum computing, [12]. While to our best knowledge it is not known at present how this type of information processing could be utilized in IP in a direct way, it is interesting to point out that some methods of manipulating matrices on a circuit-model quantum computer have been devised, [8].

Alternatively, a quantum system may be described statistically. For instance, the statistical approach is more appropriate when we discuss an ensemble of quantum systems. In such a case the state of the system is expressed via a mixed state (density matrix) ρ, which satisfies the requirements specified in Section 2.1, and its dynamics is described by the von Neumann equation

ρ ˙ = - i [ H , ρ ] .

It is not unreasonable to reflect if the von Neumann equation could be used with success to effect desirable IP operations. However, numerical evidence suggests that for a typical Hamiltonian this evolution induces strong undesirable artifacts in the image I = rho2image ( ρ , α , β ) . Of course, it remains a possibility that a useful Hamiltonian could be found in the future.

Another fundamental concept is that of quantum measurement. Quantum measurement is a physical process which aims to establish the temporal value of an observable. Observables are self-adjoint operators on . In the canonical quantum theory, a measurement corresponding to an observable, say, 𝒜 will collapse the quantum system onto a new state which must be one of the eigenvectors of 𝒜 , at the same time its corresponding eigenvalue will be the outcome of the measurement. In this setting the outcome of a measurement is again modelled statistically – it is not certain which eigenvector the system is going to collapse onto and what will be the returned value of the measurement of the observable. The expected value, however, is given by the theory precisely. When the system is characterized by a pure state | ψ , the expected value of the observable is 𝒜 = ψ | 𝒜 ψ , and when it is characterized by the mixed state ρ, it is

𝒜 = Tr ( ρ 𝒜 ) .

The right-hand side of the above formula may be viewed as a bilinear form (with arguments ρ and 𝒜 ) which prescribes a duality between states and observables. This duality gives rise to an alternative description of quantum system dynamics, which is known as the Heisenberg picture and is fully equivalent to the one outlined above (the Schrödinger picture). In the Heisenberg picture the system state, be it v or ρ remains time independent, but the observables evolve in time. Thus, the evolution of the observable 𝒜 is described by the Heisenberg equation, namely

𝒜 ˙ = i [ H , 𝒜 ] .

Again, one might consider this equation as a potential approach to IP. Indeed, the evolution would be applied to operators F 1 and F 2 where [ F 1 , F 2 ] = image2observables ( I ) and subsequently the evolved pair of observables would be back-converted to an image. However, also in this case it is not known to us whether or not a useful IP task could be accomplished in this way.

3.2 Quantum dynamical semigroups: A useful IP tool!

Modern quantum mechanics goes beyond the rudimentary framework of Section 3.1 and strives to describe the measurement process with more precision. In the prevailing contemporary approach the measurement is modeled as an interaction of two quantum systems: the one being observed and the one used to observe, i.e. the measuring apparatus. A more detailed discussion of this model may be found in [13], while here we outline the basics. Since the focus is on the observed system, the evolution is described via its state variable, the density matrix ρ. (Since in this case the quantum system is part of a larger quantum system, its state cannot be known as a pure state vector. This is rather the result of the mathematical reality of the quantum mechanical description of the system dynamics which, it is believed, reflects the physical reality.) The dynamics is now only partly driven by the Hamiltonian H of the observed system. Another contribution to the dynamics comes from the influence of the measuring apparatus, which is expressed via one or more operators, the so-called Lindblads, say, L j , j = 1 , 2 , , J . The operators L j may but need not be self-adjoint. The evolution of the state variable ρ is then given by the so-called Markovian master equation in Lindblad form, namely:

(7) ρ ˙ = [ ρ ] := - i [ H , ρ ] + j { L j ρ L j - 1 2 L j L j ρ - 1 2 ρ L j L j } .

There is also a master equation proper for the Heisenberg picture, for details see [1], Chapter 8, namely:

(8) 𝒜 ˙ = i [ H , 𝒜 ] + j { L j 𝒜 L j - 1 2 L j L j 𝒜 - 1 2 𝒜 L j L j } .

A simple yet crucial mathematical property of the master equations is that a density matrix evolving according to (7) will satisfy all the requirements for a density matrix at all times and, similarly, an observable evolving according to (8) will remain an observable at all times. A more sophisticated property of the resulting quantum semigroups is their complete positivity, [1], which we do not discuss here. Numerical simulations reveal a very complex interplay between the Hamiltonian and the Lindblads; in particular, either part can be the dominant factor in determining the evolution. We have observed that evolution of type (7) with carefully designed constituent operators and applied to ρ from [ ρ , α , β ] = image2rho ( I , α ) is indeed useful in the IP applications – see Section 4.1 for details.

3.3 Quantum state diffusion: Might it be a useful IP tool?

The derivation of the Markovian master equation in Lindblad form, (7), relies upon the theory of open quantum systems, see e.g. [1]. In the last two decades of the 20th century, researchers (Nicolas Gisin, Ian Percival, and others – see [13] for the history of the subject) developed a reinterpretation of the master equation on grounds of stochastic calculus. The resulting theory, known as the theory of Quantum State Diffusion (QSD), is extremely intriguing in that it may explain the true microscopic mechanism behind the master equation governed dynamics, [13]. It is interesting to examine if there may be advantages to applying the QSD toward improved image processing. Here, we give a brief relation of the theory. Recall that the mathematical model of Brownian motion is the so-called Wiener process, W ( t ) . Its infinitesimal increment d W ( t ) has mean zero, i.e. E [ d W ( t ) ] = 0 , and variance E [ d W ( t ) 2 ] = d t . For two independent Wiener processes d W 1 ( t ) and d W 2 ( t ) , we define d ξ ( t ) = ( d W 1 ( t ) + i d W 2 ( t ) ) / 2 . One may regard d ξ ( t ) as the Brownian motion in the complex plane. It satisfies

E [ d ξ ] = 0 , E [ d ξ 2 ] = 0 , E [ | d ξ | 2 ] = d t .

Next, let | ψ denote a pure quantum state of a quantum system. As usual, the expectation of an operator L (be it self-adjoint or not) is defined as L := ψ | L | ψ . The following observation is the cornerstone of QSD. If | ψ evolves according to the stochastic differential equation

(9) d | ψ = - i H | ψ d t + j ( L j L j - 1 2 L j L j - 1 2 L j L j ) | ψ d t + j ( L j - L j ) | ψ d ξ ( t ) ,


(10) ρ = E [ | ψ ψ | ]

satisfies the master equation (7). Note that the stochastic differential equation (9) is nonlinear, as indeed terms L j , L j depend on | ψ ( t ) . It is easily seen that the evolution preserves the normalization of ψ, i.e. ψ ( t ) = 1 for all times t. Note also that the interpretation of ρ as a mean of a plethora of pure states, (10), is characteristic for the QSD theory.

In principle (9) could be used as a framework for IP applications. Indeed, one can find ρ as usual via [ ρ , α , β ] = image2rho ( I , α ) and then find an initial ensemble of vectors | ψ to satisfy (10). Subsequently, the ensemble of state vectors can be evolved stochastically a sufficient number of times to give a good approximation of the expectation (10) at a later time, which can then be back-converted into a luminance function. While it remains to be seen if there may be any advantage to using this framework rather than the direct approach furnished by (7) in digital simulations – see additional remarks in Section 4.4 –, the stochastic interpretation is quite intriguing as a model for the physical process. It should be considered as fundamental if the algorithms discussed here find a hard-wired implementation via an engineered quantum process.

3.4 Classical simulation of quantum dynamics

From now on we will consider the model (7) or (8) with some additional assumptions. First, we will only use one Lindblad, which we will denote L. Second, we will admit the possibility that the Hamiltonian H is in fact time-dependent, and for times t [ 0 , T ] has the form

(11) H ( t ) = ( 1 - t T ) H i + t T H f ,

where H i (resp. H f ) is the initial (resp. final) value of the Hamiltonian. This model approximates the process of quantum annealing which is of particular interest in the context of adiabatic quantum computing. Since quantum computers are in a sense universal quantum systems – as, indeed, their primary purpose is to simulate arbitrary quantum processes –, it is interesting to immerse the quantum approach to image processing in this context. Of course, the model based on (11) is an idealization that does not aspire to capture the precise characteristic of quantum computing machines whose function depends on the specific hardware implementation, the nature of the machine’s interaction with the environment, the temperature, the true H ( t ) as a function of time, etc.

Note that in re-scaled time s = t / T the evolution (7) assumes the form

(12) d d s ρ = 1 [ ρ ] + s 2 [ ρ ]


1 [ ρ ] := - i T [ H i , ρ ] + T ( L ρ L - 1 2 L L ρ - 1 2 ρ L L ) ,
(13) 2 [ ρ ] := - i T [ H f - H i , ρ ] .

It is demonstrated in [17] that solutions of (12) may be found in the Taylor series form

(14) ρ ( s ) = ρ 0 + s ρ 1 + s 2 ρ 2 + s 3 ρ 3 +

with matrix coefficients

(15) ρ n = 1 n ! Q n ( 1 , 2 ) ρ 0 ,

where Q n = Q n ( 1 , 2 ) is a polynomial in variables 1 and 2 defined by recursion:[2]


Q 0 = 𝕀 ,
Q 1 = 1 ,
Q n + 1 = 1 Q n + n 2 Q n - 1 for n 1 .

As it turns out, the Taylor series converges absolutely in the entre complex plane (a fortiori for all s [ 0 , 1 ] ), and so it defines an operator-valued entire function. The recurrence formula is also the foundation of efficient numerical schemas for solving equation (12), see [17] for details. These schemas have been adapted to effect image processing as postulated in Section 3.2.

We note that this approach extends to the dual equation (8). Indeed, representing (8) in the format

(17) d d s 𝒜 = 1 [ 𝒜 ] + s 2 [ 𝒜 ] ,


1 [ 𝒜 ] := i T [ H i , 𝒜 ] + T ( L 𝒜 L - 1 2 L L 𝒜 - 1 2 𝒜 L L ) ,
(18) 2 [ 𝒜 ] := i T [ H f - H i , 𝒜 ] ,

it is easily seen that solutions are of the form

𝒜 = 𝒜 ( s ) = n s n n ! Q n ( 1 , 2 ) 𝒜 0 .

In this manner the dual equation may also be applied as a means of processing images that are represented as quantum observables.

Remark 3.1

Note that the equation solvers outlined above apply even if the Hamiltonians H i and H f are both zero. Such a “pure Lindblad” evolution is useful in some image processing tasks. Physically, it could be interpreted as evolution in which interaction with the environment overpowers the dynamics of the quantum system so strongly that all other factors, e.g. its internal structure, may be disregarded. In a complementary view the evolution with zero Hamiltonians may be interpreted as a rotating frame representation of dynamics that is in fact driven by a strongly dominating Hamiltonian. In the latter case the task of transcribing the data into the stationary frame is carried out by unitary operations.

Remark 3.2

D-Wave quantum annealers have been designed to physically implement an adiabatic evolution of the state vector. In this architecture a skillful choice of the Hamiltonians H i and H f facilitates (quantum) computation via an adiabatic process, see e.g. [6]. For this type of computation to be successful the computer must work in a regime for which the so-called Adiabatic Theorem remains valid. In particular, the quantum computing circuit must be well isolated from its environment, which in our model corresponds to a negligible Lindblad operator. There is an ongoing debate as regards the extent to which these devices implement the theorized type of evolution, see e.g. [3, 17].

Remark 3.3

It is helpful to make a clear distinction between the term quantum annealers as it is used above and a body of work extant in the Mathematical Physics literature that refers to deep mathematical analogies, first observed in [7], between the simulated annealing algorithm and certain types of quantum dynamics. Therein the point is that there exists a type of process, expressed in the formalism of von Neumann algebras, that encompasses the classical annealing as well as certain quantum processes as its special cases. This is distinct from the principle of adiabatic quantum computing mentioned above.

4 Numerical simulations

4.1 Post-processing a noisy real life image

We have applied a numerical algorithm based on equation (12) and Ansatz (14), as outlined in the previous section, to enhance a natural image, see Figures 1 and 2. By way of experimentation we have observed that the Hamiltonians H i and H f typically play an undesirable role in the process and ought to be small or completely suppressed. On the other hand, the Lindblad operator is crucial. In the experiment at hand the Lindblad was a diagonal operator with the eigenvalues being the scaled consecutive integers. However, that does not mean that the image cannot be processed in a quantum environment where the Hamiltonian dominates the Lindblad operator, see Section 3.4, Remark 3.1.

Figure 1 A noisy gray scale image (b) is enhanced by the use of the quantum flow algorithm with T=18${T=18}$ (c) and T=36${T=36}$ (d). The clean original image (a) comes from a photograph of the Haymarket Memorial Clock Tower in Leicester, England – a view at the statue of a 15th/16th century merchant William Wyggeston.)

Figure 1

A noisy gray scale image (b) is enhanced by the use of the quantum flow algorithm with T = 18 (c) and T = 36 (d). The clean original image (a) comes from a photograph of the Haymarket Memorial Clock Tower in Leicester, England – a view at the statue of a 15th/16th century merchant William Wyggeston.)

Figure 2 A noisy gray scale image of an architectural ornament is enhanced by the use of the quantum flow algorithm, also the clean image has been subdued to the same process. Note that (with the particular setting of the parameters) the process has not introduced any drastic distortions.

Figure 2

A noisy gray scale image of an architectural ornament is enhanced by the use of the quantum flow algorithm, also the clean image has been subdued to the same process. Note that (with the particular setting of the parameters) the process has not introduced any drastic distortions.

4.2 Pure Lindblad flow as a filter

While quantum dynamics is typically counterintuitive, in the particular case highlighted in Section 4.1 it is possible to analyse its effects explicitly and unequivocally. Indeed, let the Lindblad operator be given as

L = ( l 1 0 0 0 l 2 0 0 0 l 3 ) ,

and let ρ = [ ρ m , n ( t ) ] evolve according to the pure-Lindblad evolution equation

(19) ρ ˙ ( t ) = [ ρ ] := L ρ L - 1 2 L L ρ - 1 2 ρ L L .

An explicit matrix calculation shows that the solution is given by

ρ m , n ( t ) = ρ m , n ( 0 ) exp ( l m l ¯ n - 1 2 | l m | 2 - 1 2 | l n | 2 ) t .

Thus, in the particular case l n = c n , where c is a real normalizing constant, the Lindblad evolution executes a filtering operation:

(20) ρ m , n ( t ) = ρ m , n ( 0 ) exp [ - c 2 2 ( m - n ) 2 t ] .

Recall that ρ encodes the Fourier coefficients of the image. Although a filter of this type may well have been invented outside the scope of quantum dynamics, we are not aware of it being reported in image processing literature.

4.3 Entropy production

It is well known, [4], that the evolution of an open quantum system given by (7), a fortiori (19), results in entropy production. Namely, let τ be an arbitrary stationary state, i.e. a state that satisfies [ τ ] = 0 . Then the relative entropy defined by

S ( ρ ( t ) τ ) = Tr ( ρ ln ρ ) - Tr ( ρ ln τ )

is a non-increasing function of time, i.e. the entropy production rate

σ ( t ) = - d d t S ( ρ ( t ) τ )

is nonnegative for all times t > 0 . It is easily seen, for details see e.g. [4], that σ ( t ) is determined via a solution ρ = ρ ( t ) of (7) as follows:

(21) σ ( t ) = - Tr ( [ ρ ] ln ρ ) + Tr ( [ ρ ] ln τ ) .

It is interesting to characterize σ ( t ) in the context of dynamics given by (19). Below we calculate the special case of (20) with 2 × 2 matrices and show that the entropy production rate decays exponentially in time. It means that most of the entropy dispensed in the entire process is in fact dispensed very rapidly at the beginning. Of course such a result is not unexpected in light of the exponential decay of the off-diagonal terms of ρ ( t ) . In consequence it is natural to conjecture exponential decay of the entropy production rate for the dynamics of type (19) in any number of dimensions. This leads one to believe that image filtering via such processes (with suitable Lindblads) is time efficient.

Calculation that follows is facilitated by the Pauli matrices:

σ x = ( 0 1 1 0 ) , σ y = ( 0 i - i 0 ) , σ z = ( 1 0 0 - 1 ) .

Namely, let the initial value of the density matrix be

ρ ( 0 ) = ( a x + i y x - i y b ) = 1 2 𝕀 + x σ x + y σ y + a - b 2 σ z ,

where a , b > 0 , a + b = 1 . Next, let us assume for simplicity c = 2 and note that (20) implies

(22) ρ ( t ) = 1 2 𝕀 + e - t x σ x + e - t y σ y + a - b 2 σ z = : 1 2 ( 𝕀 + T ) .

An elementary calculation shows that for ρ ( t ) to admit two positive eigenvalues whose sum is 1 it is necessary and sufficient that Δ = 4 ( x 2 + y 2 ) e - 2 t + ( a - b ) 2 < 1 for all times t, which is equivalent to the condition x 2 + y 2 < a b . Note that in the present notation

(23) [ ρ ] = d d t ρ ( t ) = - e - t x σ x - e - t y σ y .

Furthermore, note that (22) implies

(24) ln ( 2 ρ ( t ) ) = ln ( 𝕀 + T ) = k = 1 ( - 1 ) k + 1 k T k ,

and, by a straightforward calculation,

T k = { Δ k / 2 𝕀 , k even , Δ ( k - 1 ) / 2 T , k odd .

Next, we observe that the series in (24) converges for Δ < 1 , and

(25) ln ( 2 ρ ( t ) ) = 1 2 ln ( 1 - Δ ) 𝕀 + 1 2 Δ ln 1 + Δ 1 - Δ T .

In order to evaluate the entropy production rate, we take τ = ( 1 / 2 ) 𝕀 to be the reference density. Clearly, we have [ τ ] = 0 . We also note that Tr [ ρ ] = 0 , so that Tr { [ ρ ] ln τ } = 0 . Since in addition ln { 2 ρ } = ln ( 2 𝕀 ) + ln ρ , formula (21) assumes a simplified form

(26) σ ( t ) = - Tr ( [ ρ ] ln ρ ) = - Tr ( [ ρ ] ln ( 2 ρ ) ) .

Now, combining (23), (25) and (26), we readily obtain

σ ( t ) = 2 e - 2 t ( x 2 + y 2 ) 1 Δ ln 1 + Δ 1 - Δ 2 e - 2 t ( x 2 + y 2 ) 1 | a - b | ln 1 + | a - b | 1 - | a - b | ,

where the right-hand side describes the asymptotic behavior. In summary, the entropy production rate decays exponentially.

4.4 Conditions adverse to image enhancement

We have tested the effect of quantum evolution on the noisy image under a variety of conditions. While a longer more systematic study is needed to fully characterize the set of conditions conducive as well as conditions adverse from the viewpoint of image analysis, we give here a few preliminary remarks:

  • We have simulated the quantum flow in the QSD setting, i.e. equation (9). We find that this approach is likely not suitable as a basis for numerical image processing algorithms. Indeed, in light of the fact that the density matrix ρ is given as the expectation of random projectors it is necessary to generate a plethora of solutions to the stochastic differential equation (9)). This turns out to be prohibitively costly computationally so that, in the end, we have not seen sufficient convergence. However, we haste to add, the QSD schema remains attractive as a model for real world quantum state dissipation.

  • We have tested the effect of adjusting the strength of the Hamiltonian (11) be it stationary (e.g. H f = 0 ) or evolving. We have found that when the Lindblad operator L is diagonal, as in Section 4.2, strong Hamiltonians will introduce unacceptable artifacts into the image. However, it is not necessary that the Hamiltonians be zero for suppression of these artifacts. Rather, it suffices that the Lindblad dominates the Hamiltonians in their characteristic magnitudes.

  • We have tested quantum-dynamics based image evolution in the case when the image is represented via a pair of observables (6). (Of course, in this instance the numerical experiment requires discretization of the variables F 1 and F 2 .) Unfortunately, this type of processing introduces an undesirable diagonal “seam” which would have to be dealt with in a post-processing stage. Again, future research may indicate the means for removing this type of an artifact. However, the main issue is that a hard wired quantum process, if it is successfully implemented, will be free of the numerical artifacts.

5 Conclusions

We have examined the possibility of effecting a desirable natural image transformation by means of simulated quantum dynamics. Our findings suggest that functionalities such as denoising may be effected in quantum environments in conditions of decoherence. In brief, one established scenario for denoising is that of a strong interaction of a pre-designed environment with the quantum system that holds the image information encoded in a quantum variable. In particular, if a quantum device were engineered to effect such results, its regime of operation would not be conducive to standard quantum computing as the latter requires low decoherence.

One of the philosophical puzzles of modern physics is the emergence of geometry in a principally quantum-mechanical world. Although literature on the subject is limited, some intriguing mathematical speculations have been proposed, see e.g. [2]. It is interesting to view the technique discussed here as a contribution toward that theme, particularly as regards the role of decoherence in stabilizing classical (e.g. geometric) information. Finally, the technique of representing spatial data as quantum observables is further developed in the recent work [16].

A Optimal selection of parameters in the construction of a density matrix

How can one estimate the constant α used in (2) that will ensure positive semi-definiteness of ρ? It follows from the theorem of Lidskii, see [9, Chapter II.5], that the eigenvalues of ρ are shifted with respect to the eigenvalues of X + X by a vector, say, V that lies in the convex hull of the set of vectors obtained by permuting the components of [ α D ˇ 1 , , α D ˇ N ] . Let us examine the consequences of this fact.

First, let

I p = ( k , l | I k , l | p ) 1 / p

be the p norm of the luminance function I 0 . We will only make use of this norm in the two cases p = 1 , 2 . We also recall the well-known fact that the induced norm of the matrix operator I (denoted I ) is bounded above by the 2-norm (known as the Hilbert–Schmidt norm), i.e. I I 2 . Using this, it is easy to estimate the eigenvalues of X + X . Indeed, X + X is self-adjoint and so diagonalizable. Let λ be an eigenvalue. We have

| λ | X + X X + X 2 I ^ 2 = I 2 ,

which shows that all eigenvalues of X + X lie in the interval [ - I 2 , I 2 ]

Next, we take a closer look at D ˇ . It is easily seen that the following identity holds (best communicated in the pseudo-code notation):

D ˇ N - j = sum(diag( rot 180 I , - j ) ) + sum(diag( rot 180 I , N - j ) ) for j = 0 , 1 , , N - 1 .

In other words, the consecutive terms of D ˇ are obtained as follows: First, stack two copies of rot 180 I one on top of another, thus forming a rectangular 2 N × N matrix. Next, compute sums of diagonal lines starting at the upper left corner ( 1 , 1 ) , then repeat for the diagonal starting at ( 2 , 1 ) , then the diagonal starting at ( 3 , 1 ) , etc.[3] In particular,

k D ˇ k = m n I m , n .

Hence, one can expect D ˇ k I 1 / N .

In light of these observation, in order to ensure that ρ is positive semi-definite, it suffices to take α such that each component of [ α D ˇ 1 , , α D ˇ N ] is at least as large as I 2 . Indeed, this will ensure that the shift by V will move the entire interval [ - I 2 , I 2 ] inside [ 0 , ) . It means that typically it should be sufficient to take

α N I 2 I 1 .

Note that by the Cauchy–Schwarz inequality I 2 / I 1 1 / N , and the ratio of the two-norm to the one-norm is closer to one for flatter (closer to constant) functions I. If the luminance function is represented via floating-point numbers in the interval [ 0 , 1 ] , then trivially I 2 2 I 1 , and so α N / I 2 should suffice.


I wish to express immense gratitude to Dr. Paul Babyn – my interest in Image Processing would not have been renewed without his help, support, and encouragement. While Paul felt his contributions did not warrant the status of a co-author, they have been numerous and important. From the very moment of inception he encouraged me to develop the ideas now expounded in this article. He also helped evaluate the quality of numerical results, provided guidance on modern IP literature, and let me lean on his expertise in radiology as I was trying to understand the scope of applicability for IP results. I am equally grateful to Dr. Alexandre Zagoskin for stimulating conversations, enlightening comments on Quantum Engineering, as well as encouragement and support for this pursuit. I am also indebted to Dr. Javad Alirezaie, who helped me better understand the state of the art image processing, Dr. Martin Laforest for many interesting comments on quantum information processing, and to Dr. Gordon Sarty for helpful and encouraging feedback on this work. A preliminary report of this work was given at the International Focus Workshop on Metamaterials and Quantum Criticality, Loughborough, Leicestershire, UK, March 25–26, 2014, and at the 11th World Congress of Computational Mechanics and 5th European Conference on Computational Mechanics, Barcelona, Spain, July 20–25, 2014. I am also grateful to the referee for pointing out additional references.


[1] Alicki R. and Fannes M., Quantum Dynamical Systems, Oxford University Press, Oxford, 2009. Search in Google Scholar

[2] Baird P., Emergence of geometry in a combinatorial universe, J. Geom. Phys. 74 (2013), 185–195. Search in Google Scholar

[3] Boixo S., Rønnow T. F., Isakov S. V., Wang Z., Wecker D., Lidar D. A., Martinis J. M. and Troyer M., Evidence for quantum annealing with more than one hundred qubits, Nat. Phys. 10 (2014), 218–224. Search in Google Scholar

[4] Breuer H.-P. and Petruccione F., The Theory of Open Quantum Systems, Oxford University Press, Oxford, 2002. Search in Google Scholar

[5] Chan T. F. and Shen J., Image Processing and Analysis: Variational, PDE, Wavelet, and Stochastic Methods, SIAM, Philadelphia, 2005. Search in Google Scholar

[6] Farhi E., Goldstone J., S. , Gutmann , Lapan J., Lundgren A. and Preda D., A quantum adiabatic evolution algorithm applied to random instances of an NP-complete problem, Science 292 (2001), 472–476. Search in Google Scholar

[7] Frigerio A., Simulated anenaling and quantum detailed balance, J. Stat. Phys. 58 (1990), 325–354. Search in Google Scholar

[8] Harrow A. W., Hassidim A. and Lloyd S., Quantum algorithm for linear systems of equations, Phys. Rev. Lett. 103 (2009), Article ID 150502. Search in Google Scholar

[9] Kato T., Perturbation Theory for Linear Operators, Springer, New York, 1966. Search in Google Scholar

[10] Krzakala F., Mézard M., Sausset F., Sun Y. F. and Zdebrová L., Statistical physics based reconstruction in compressed sensing, Phys. Rev. X 2 (2012), Article ID 021005. Search in Google Scholar

[11] Marsousi M., Abhari K., Babyn P. and Alirezaie J., An adaptive approach to learn overcomplete dictionaries with efficient numbers of elements, IEEE Trans. Signal Process. 62 (2014), 3272–3283. Search in Google Scholar

[12] Nielsen M. A. and Chuang I. L., Quantum Computation and Quantum Communication, Cambridge University Press, Cambridge, 2000. Search in Google Scholar

[13] Percival I., Quantum State Diffusion, Cambridge University Press, Cambridge, 1998. Search in Google Scholar

[14] Sowa A., Factorizing matrices by Dirichlet multiplication, Linear Algebra Appl. 438 (2013), 2385–2393. Search in Google Scholar

[15] Sowa A., The Dirichlet ring and unconditional bases in L 2 [ 0 , 2 π ] , Funct. Anal. Appl. 47 (2013), 227–232. Search in Google Scholar

[16] Sowa A., Encoding spatial data into quantum observables, preprint 2016, Search in Google Scholar

[17] Sowa A. P., Everitte M. J., Samson J. H., Savilev S., Zagoskin A. M., Heidel S. and Zúñiga-Anaya J. C., Recursive simulation of quantum annealing, J. Phys. A 48 (2015), Article ID 415301. Search in Google Scholar

Received: 2016-2-15
Accepted: 2017-1-10
Published Online: 2017-1-28
Published in Print: 2017-3-1

© 2017 by De Gruyter