## Archive for the ‘**Maths**’ Category

## How does an absorbing boundary modify the propagator of a stochastic process?

A recurring theme throughout my blog are stochastic processes, widely used to model quantities under the influence of both deterministic evolution laws and random fluctuations. A complete solution for a stochastic process is its *propagator*: The probability density for the final state at some time given an initial condition at time . The propagator depends on the boundary conditions imposed on the stochastic process. For example, an absorbing boundary is useful for describing extinction in a stochastic model of population dynamics: When the population of a given species reaches zero, it cannot be re-introduced.

Typically, however, the “free” propagator for a stochastic process — i.e. the propagator without any boundaries — is much easier to determine than the propagator with any given boundary condition.

However, common sense tells us that the free propagator for a stochastic process contains already all there is to know about it. So I have wondered for quite a while if and how one can recover e.g. the propagator with an absorbing boundary, knowing only the free propagator. I found the solution only quite recently in this preprint on arXiv, and think that it is interesting enough to explain it here briefly.

## The result

Consider a one-dimensional, stationary, Markovian stochastic process. Let us denote its *free propagator* by

This is the density for the final value at time , given an initial value at time . Stationarity of the process implies that is a function of the time difference only. Define the Laplace transform of the free propagator by

Now let us add an absorbing boundary at . Then, the propagator

is the density for the final value at time , given an initial value at time , *without touching the line between and *. **The upshot of this post is the following simple relationship between the Laplace transforms of the free propagator and of the propagator with absorbing boundary :**

(1)

In case the denominator on the right-hand side is zero, one may have to perform some limiting procedure in order to make sense of this formula. We can now check how nicely it works for simple cases, like Brownian motion with drift, where the result is known by other means.

## Example: Brownian motion with drift

Let us consider an example, the Brownian motion with drift . Its stochastic differential equation is

where is white noise with .

The free propagator is a simple Gaussian

Its Laplace transform is given by:

To compute , let us assume , . Note this is without loss of generality, since if and are on opposite sides of the absorbing boundary the propagator is zero. Applying the result (1), we have

The inverse Laplace transform of this expression is:

This is, of course, exactly the same result (the so-called **Bachelier-Levy formula**) which I had already discussed earlier in a guest blog post at Ricardo’s blog, using slightly different methods (a variant of the Cameron-Martin-Girsanov theorem). Note that some signs which were wrong in the original post are corrected here.

## The proof

Now let us go back to the case of a general stochastic process, and discuss how formula (1) can be proven. Let us assume that the stochastic process in question has a path-integral representation, i.e. that

for some action .

Now, to add an absorbing boundary at we weigh each crossing of by a factor in the path integral, and then let . This effectively retains only those trajectories in the path integral, which never touch . We thus have:

.

Now, to evaluate the path integral, we expand it in a series in (or, in other words, in the number of crossings of ). We get:

.

Note that here I cancelled the factor , coming from the Taylor expansion of , against the possibilities of ordering the times . This is why I can assume them to be in ascending order above. Also, the functions actually set all .

Taking Laplace transforms, the time integrals (which are a -fold convolution) simplify to a simple product:

.

Now, one easily obtains the claimed result:

.

This computation can, of course, be generalized to more complicated situations like several boundaries, higher dimensions, etc. For details see again the original preprint by C. Grosche. Have fun and let me know if you have questions!

## Tricks for inverting a Laplace Transform, part IV: Substitutions

After a few less technical posts recently, I now continue the series of articles on tricks for the inversion of Laplace transforms. You can find the other parts here: part I (guesses based on series expansions), part II (products and convolutions), part III, and part V (pole decomposition).

Today’s trick will be **variable substitutions**. Let’s say we need to find so that

, (1)

where , and are constants. This example is not as absurd as it may seem; it actually came up recently in my research on a variant of the ABBM model.

There is a general formula permitting one to invert the Laplace transform above, in terms of an integral in the complex plane:

.

This so-called Fourier-Mellin integral runs along a contour from to . can be chosen arbitrarily, as long as the integral converges and all singularities (poles, branch cuts, etc.) of lie on the right of it. Note that my convention for the Laplace variable is the convention used for generating functions in probability theory, which has just the opposite sign of of the “analysis” convention for Laplace transforms.

The fact that our Laplace transform is written entirely in terms of the function suggests applying a variable substitution. Instead of performing the Fourier-Mellin integral over , we will integrate over . We then have:

. (2)

Solving the definition of for , we get

and

.

Inserting this into (2), we have

.

The only singularity of the integrand is at , for . We can thus choose the contour to go parallel to the imaginary axis, from to . `Mathematica`

then knows to how to evaluate the resulting integral, giving a complicated expression in terms of Kummer’s confluent hypergeometric function . However, a much simpler expression is obtained if one introduces an auxiliary integral instead:

.

`Mathematica`

knows a simple expression for it:

.

is Tricomi’s confluent hypergeometric function, which is equivalent to Kummer’s confluent hypergeometric function but gives more compact expressions in our case.

Using this auxiliary integral, (2) can be expressed as

.

Simplifying the resulting expressions, one obtains our **final result for **:

(3)

Equivalently, the confluent hypergeometric functions can be replaced by Hermite polynomials:

For complicated Laplace transforms such as these, I find it advisable to check the final result numerically. In the figure below you see a log-linear plot of the original expression (1) for , and a numerical evaluation of , with given by (3). You can see that they match perfectly!

## Random matrix theory and the Coulomb gas

*Today I have the pleasure of presenting you a guest post by Ricardo, a good physicist friend of mine in Paris, who is working on random matrix theory. Enjoy!*

After writing a nice piece of hardcore physics to my science blog (in Portuguese, I am sorry), Alex asked me to come by and repay the favor. I am happy to write a few lines about the basis of my research in random matrices, and one of the first nice surprises we have while learning the subject.

In this post, I intent to present you some few neat tricks I learned while tinkering with Random Matrix Theory (RMT). It is a pretty vast subject, whose ramifications extend to nuclear physics, information theory, particle physics and, surely, mathematics as a whole. One of the main questions on this subject is: given a matrix whose entries are taken randomly from a known distribution, what would the distribution for its eigenvalues be?

The applications of such a question are countless. We might be interested in the energy levels of a random hamiltonian, the information transmitted through a random noisy channel, the transmission probability of an electron through an open quantum system; every eigenvalue problem, when the effects of the environment are taken into account, must answer this most fundamental question in RMT. To do so, first we restrict our analysis to symmetric (resp. hermitian) matrices. We have two reasons for it: firstly, most applications in physics involve this kind of matrix (hamiltonians, S-matrices), secondly, they are much easier to analyze, since they can always be diagonalized by a orthogonal (resp. unitary) matrix.

To enable a complete analysis, I must introduce an extra assumption: our matrix must have a probability density function (p.d.f) invariant with respect to orthogonal (resp. unitary) conjugation. This is not such a narrow class of matrices as one might think, we can construct this kind of matrix with nearly any infinitely divisible probability. Roughly, we want to be able to write:

where is a diagonal matrix.

However, we should not get too excited with this property, as to write right away. We are always talking about a probability **density** function, the computation of the actual probability requires an integral over a certain space, and the passage from the matrix to its eigenvalues requires a variable change, hence we must calculate the Jacobian of such a transformation.

Let us take the simplest and perhaps most useful example: the Gaussian real (resp. complex) symmetric (resp. hermitian) matrix whose entries are i.i.d. normal variables with zero mean and variance 1. We can verify that this matrix is invariant with respect to orthogonal conjugation (actually, we can prove that the Gaussian matrix is the **only** matrix invariant with respect to orthogonal conjugation and whose entries are independent). Hence, we expect the eigenvalue probability to be:

The value for will be discussed soon, is a normalization constant. I will not derive the Jacobian term and I will refer to [1] for the complete demonstration, not only of this expression, but also of the uniqueness theorem mentioned above. Even though the entries of our matrix are independent, the p.d.f. of its eigenvalues its highly coupled due to the Jacobian:

is one of three possible values: if our initial matrix is real, . If it is a complex matrix, . If is a quaternion matrix, . This is called **Dyson’s threefold way**.

Since we have the eigenvalues p.d.f., it should be easy to calculate all quantities involved in a eigenvalue problem. We could be interested in the probability of having an eigenvalue larger than a threshold [2], or the chance of not finding any eigenvalue at all inside an interval. The calculation necessary would be simply the -dimensional integral of the p.d.f. multiplied by a function, but, oh dear, what a probability function!

The Jacobian will make most variables changes futile, all attempts of uncoupling the equations are faced with the product . One standard approach would be to interpret the Jacobian as a Vandermonde determinant (up to a sign), rearrange rows to compose polynomials and adjust them as to create orthogonal functions, whose integral would preserve only the “diagonal” terms of our calculations. I will not explain this technique in detail, even though it is quite beautiful. I will travel another path.

Let us write the whole probability as an exponential. We already have the Gaussian term as an exponential, let us raise the Jacobian to the exponent in the most standard way, we write:

where . If we write , the analogy is complete, we can treat this probability density function as a system analyzed by the Boltzmann Canonical Ensemble! Hence, we import all our techniques and tools from statistical mechanics to solve this ugly, nasty piece of probability density function.

It is interesting to remark that this Hamiltonian has a natural physical interpretation. If you recall your electrodynamics course, you might find, deep in your memory, the derivation of the Coulomb potential for two point charges in 2D: if we say their position is and , we have . Our Hamiltonian is the description of a 2D gas under a quadratic potential, or, as we call it, **Dyson’s Coulomb gas**. We can actually infer the equilibrium position of these electrons, if we consider the quadratic potential and the mutual logarithmic repulsion:

It is no surprise that, for the large limit, the p.d.f. of the electrons, and hence that of the eigenvalues, is the Wigner Semicircle Law. To deduce it, the next natural steps are (I am not going to detail them any further): replace by its density: , write the Hamiltonian in terms of the density, replace the integral over the eigenvalues by a functional integral over the possible densities , use a saddle-point approximation to argue that only the “best” density will matter for our calculations, derive functionally the Hamiltonian and solve a Tricomi equation to obtain the optimal value for , the equilibrium position of the electrons. Surely, we suppose that is large for the saddle-point approximation, results for small require a different approach.

To wrap it up, we can avoid treating a highly coupled integral by sending the nasty coupled part to the exponent, interpret the system as a Canonical Ensemble, and treat it in the thermodynamical limit to find the density of eigenvalues. This beautiful treatment was popularized by Dyson and revisited many times in the last decade. RMT is a very active subject, for physicists and mathematicians, and I hope that, with these few lines, you have been able to see a little but of what I see on it, and to appreciate, for a change, a place where physics has come to help mathematics on a hard day of work.

[1] M. L. Mehta, *Random Matrices*, Third Edition, Elsevier, 2004.

[2] S. N. Majumdar, C. Nadal, A. Scardicchio, P. Vivo, *How many eigenvalues of a Gaussian random matrix are positive?*, Phys. Rev. E 83, 041105 (2011)

## A quick introduction to the Martin-Siggia-Rose formalism

**The Martin-Siggia-Rose (MSR) formalism is a method to write stochastic differential equations (SDEs) as a field theory formulated using path integrals.** These are familiar from many branches of physics, and are simple to treat perturbatively (or, in some cases, even to solve exactly).

We start from an SDE of the form:

,

where is a deterministic function and is Gaussian noise with an arbitrary correlation function :

.

MSR recognized that observables, averaged over solutions of this SDE, can be written as a path integral with a certain action:

,

where

.

A useful **special case** is an automonous (time-independent) SDE with **multiplicative noise**:

,

where is Gaussian white noise. The resulting action for the path integral is local in time:

.

I will now show a few concrete examples, and then explain how this mapping is derived.

## Example: Ornstein-Uhlenbeck process

The Ornstein-Uhlenbeck process is one of the few stochastic processes treatable analytically. Its SDE is:

,

where (else the process blows up) and is white noise:

.

In this section, I will show that the **generating functional** for this Ornstein-Uhlenbeck process is given by:

.

To show this, we apply the MSR formula to express the generating functional as a path integral:

.

This path integral can be solved exactly. The key observation is that the exponent is linear in , and hence one can integrate over this field. Using , this gives (up to a normalization factor ) a -functional:

.

The factor remains since it is the only term in the action independent of .

The -functional fixes in terms of :

.

From this we get

The normalization factor is seen to be by imposing that for we must have .

Plugging this in gives the generating functional, as claimed above:

.

One simple consequence is the **exponential decay of the connected correlation function**:

.

But of course, the full generating functional tells us much more – in some sense, it is the exact solution of the SDE!

*Note added 01/12/2012: Of course, the generating functional can be written down immediately once one knows the two-point correlation function and the fact that the Ornstein-Uhlenbeck process is Gaussian. However, the latter fact is not immediately obvious (and was not used in our derivation), in fact one can see our MSR calculation as one way to prove Gaussianity for the OU process.*

## Derivation of the Martin-Siggia-Rose formula

Now that we’ve seen what the MSR path integral formulation is useful for, let us see how it can be derived. The average over an observable can be written as a path integral with the SDE enforced through a -functional:

.

Now, we rewrite the -functional using the formula . Since the -functional is a product of a -function at all points in space, effectively we are introducing an auxiliary field :

.

Now, the only thing that depends on the noise is the factor . Since is Gaussian, the average of this factor can be simply expressed as

Altogether we obtain

which is just the MSR formula claimed above.

Amazingly, I did not yet find a clear derivation of this path integral representation on Wikipedia or anywhere else on the web… So I hope this will be useful for someone starting to work with the formalism!

If there is interest, I may write another blog post in the future with more examples and explanations e.g. on the physical role of the auxiliary field .

## Tricks for inverting a Laplace Transform, part III

*This is a continuation of the series of articles on Laplace transforms… You can also have a look at part I (guesses based on series expansions), part II (products and convolutions), part IV (substitutions), and part V (pole decomposition).*

This time we will cover the following inverse Laplace problems:

The basic result that can be used for such problems is the Laplace transform of the Bessel function:

Now, using

,

and by rescaling , one finds that is solved by:

Have fun!

## An interesting integral equation

Recently, I discussed an interesting problem with a friend of mine. Basically, the objective was to find an approximation for the solution of a linear equation system in the limit of .

It turned out that in this limit the solution can be obtained by solving the following integral equation

,

where is known and is to be found. Now, such integral equations are notoriously difficult to solve (in particular, much more difficult than differential equations). An excellent resource (which contains most known exact solutions) is the EqWorld page on integral equations (and the “Handbook of Integral Equations” by EqWorld’s authors, which you can find by googling). Unfortunately, it does not help in this case (the closest thing one finds there are the eigenfunctions for the integral operator , which, however, turn to to be competely unrelated to those of ).

One angle of attack to our equation is the observation that plugging in a polynomial gives a polynomial (and vice versa). The mapping of coefficients is linear, and by inverting its matrix one can explicitely write the solution for any given polynomial . Looking at the resulting coefficients polynomials of small degree one can guess the following expression for the solution (which is then easily verified to hold for general holomorphic ):

Here, the contour has to be chosen in such a way that all singularities of the integrand are inside . For example, one can choose a circle with radius . For non-holomorphic , this formula still works but the choice of integration contour is more complicated.

In the problem which motivated that equation, we had given for and wanted to obtain for . In this case, the contour integral formula can be rewritten as

This integral is assumped to be regularized by taking the Cauchy principal value. Computing this (e.g. numerically) now only requires for , so it is no longer necessary to assume that can be extended analytically to . A similar expression can also be found for where .

If somebody finds another use or application for this pretty result, please drop me a comment below 😉