Wednesday, June 24, 2009

CS: Testing the Non-Negative Nullspace Property for the Linear Transport Equation ?

The linear transport equation or the Linear Boltzmann Equation is hard to solve even in the linear case. You may recall my mention of a presentation by Hongkai Zhao on A fast solver for radiative transport equation and applications in optical imaging where I was not sure that the RIP was such a good benchmark. Ever since I have wondered how other conditions recently presented might be a little more appropriate to show that solving the transport equation can be performed using Linear Programming. Here is the idea. Let us say we are trying to solve the following problem:

H f = S (1)

where H is a linear transport operator, S is the known source function and the unknown is function f. The reason this type of problem is not readily solved is that the dimension of the problem is large and therefore any type of sampling/discretization yield very large matrices. In production codes, when one solves the real problem, one never produces the discretized version of H, rather the element of the matrix are computed on the fly. Back to our problem, one of the first thing to solve for is the attendant homogenous equation

H f = 0

There have been a number of publications on the subject yielding families of functions \psi with index i (see [1]) such that:

H \psi_i = 0

\psi_i represents elements of the nullspace of operator H or simply the solution of the homogenous equation (1).

The physics of the problem requires the solution f to (1) to be sparsely decomposed in a new basis \phi_i made up of positive functions such as splines. One can also expect the solution for (1) to involve a set of positive coefficients when projected on the new basis \phi_i (this also probably require the use of a multiscale set of splines not unlike the ones studied by Thierry Blu and Michael Unser)

Sparsity and Positivity, mmhhh...we have seen these elements before. More precisely, in a presentation on the nullspace property for positive signals entitled Sparse Recovery of Positive Signals with Minimal Expansion presented by Alex Dimakis, a joint work with M. Amin Khajehnejad , Weiyu Xu, Babak Hassibi (the attendant paper is Sparse Recovery of Positive Signals with Minimal Expansion ).

So with f = \Sigma w_i \phi_i = [...phi...] [x], and x sparse, we are now trying to solve for:

H [...phi...] [x] = S (2)


min||x||_0 s.t H[...phi...][x] = S

This expression is equivalent to

min||x||_0 s.t H[...phi...][x] = S


w in the nullspace of H[...phi...], w has at least k negative elements.

The nullspace condition can therefore be found when one finds the w's. They can be found to be as follows:



[...psi...]^T [...phi...][w] = I


[w] = ([...psi...]^T [...phi...])^(-1)

([...psi...]^T [...phi...]) is the matrix with element equal to \int(phi_i * psi_j).

Every column of that w matrix should then be inspected and one should set the smaller number of elements of every column as negative elements and put the larger number of elements as positive (since both \psi_i and \psi_i are elements of the nullspace of H) then one should check that "...the sum of every n − k nonnegative elements should be greater than the absolute sum of the rest...". It looks easy :-)

[1] Linear Transport Theory by Kenneth M. Case and Paul F. Zweifel

Credit: NASA, Sarychev Peak Eruption, Kuril Islands as taken by the astronauts from the International Space Station 12 days ago.

No comments: