User rcompton - MathOverflow

Fixed point theorem on graphs?

2013-02-09T22:46:09Z

Originally posted here: http://math.stackexchange.com/questions/276167/fixed-point-theorem-on-graphs

I have a graph $G=(V,E)$ where to each vertex $v$ I have associated a value, $\hat{v}$ (ie I have a "network" in the terminology here http://snap.stanford.edu/snap/index.html ).

Let $\phi : \hat{V} \rightarrow \hat{V}$ be the function which takes the value associated to each node and replaces it with median of the values of the adjacent nodes.

Empirically, iterating $\phi$ converges. Why?

edit: The graph is large and follows this model: http://en.wikipedia.org/wiki/Barab%C3%A1si%E2%80%93Albert_model . Also, for ~1% of nodes, $\phi$ is the identity. There is some related work here: http://www.cs.cmu.edu/~zhuxj/pub/CMU-CALD-02-107.pdf

edit: Looks like I can just replace the $P(j \rightarrow i)$ in the label propagation paper with $P(j \rightarrow i) = 1\; if \; \text{j is the median, }0 \text{ else}$, then I can copy their convergence result.

edit: Wait, no. If I do that then P changes each step so it's not a direct copy...

Answer by rcompton for 2D Problems Which are Easier to Solve in 3D

2012-09-09T04:44:19Z

This is essentially the idea behind level set methods (cf http://en.wikipedia.org/wiki/Level_set_method ).

There are several situations where one needs to study the behavior of a dynamic surface with complicated topology. In fire simulation, one needs to track the motion of an air/fuel interface which is often not connected. In image segmentation, one needs to move the boundary between inner/outer regions until a steady state is reached. In fluid mechanics, breaking waves detach from the main body of water.

Solving differential equations on surfaces with complicated topology is difficult numerically. It turns out that when one represents these surfaces as level sets of a higher dimensional function these connectivity problems disappear and higher quality simulations are possible.

Quantifying the amount of structure in a data set via random matrix theory

2012-08-22T23:16:35Z

Given a data matrix, $M \in \mathbb{R}^{n \times p}$, I am interested in methods quantifying the amount of structure in present in $M$.

I've found a few approaches, but I would like to learn more about what is available and what can be done given new developments in random matrix theory. I understand that the phrase "quantifying the amount of structure" is quite vague, any references that can clear that up would be appreciated.

Currently, I know about the following methods:

-- Compare the leading singular value of $M$ against the leading singular value of a shuffled version of $M$ (cf. http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.0030160).

Here, one determines whether or not $M$ has structure by comparing the leading singular value of $M$, $\sigma_{max}(M)$, with the leading singular value of a matrix, $\tilde{M}$, obtained by permuting the entries of $M$. If $| \sigma_{max}(M) - \sigma_{max}(\tilde{M}) | / \sigma_{max}(M) < 0.15$ then $M$ is said to contain no structure. Otherwise, the leading rank-1 approximation to $M$ is subtracted off and the process is repeated until no structure remains.

-- Compare the leading singular value of $M$ against the Tracy-Widom distribution (cf. http://projecteuclid.org/DPubS?service=UI&version=1.0&verb=Display&handle=euclid.aos/1009210544 and http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.0020190 )

This is identical to the previous test with the exception that we test against the Tracy-Widom distribution rather than the singular values of a shuffled version of $M$. I suppose this test is more rigorous but will only tell whether or not the entries of $M$ are Gaussian ( though I get the impression that newer results improve on this, eg http://arxiv.org/pdf/0906.0510v10.pdf ).

-- Check how well the correlation eigenvalues of $M$ agree with the semicircular law (cf. http://arxiv.org/pdf/cond-mat/0108023v1.pdf )

Here, the authors study the structure of several time series by observing how many eigenvalues lie within Wigner's semicircle. It would be nice if the number of significant principle components one gets from the methods above is equal to the number of eigenvalues which lie outside the semicircle. I have no idea if this is true.

Convolutive noise removal

2012-06-03T03:30:43Z

I have the time domain signal $$ u_o(t) = u(t)e^{-t/\tau}\eta(t) + \sigma(t) $$ where $\tau$ is known, $\eta$ is non-Gaussian noise, and $\sigma$ is Gaussian noise. The distribution of $\eta(t)$ is known, but only numerically. I also have prior knowledge that $u(t)$ is a sum of a small number of sinusoids. How can I recover $u(t)$ from $u_o(t)$?

In the case where $\eta$ is not present, I can Fourier transform to obtain: $$ \hat{u}_o(\xi) = \hat{u}(\xi)*l(\xi) + \sigma(\xi) $$ where $l$ is a Lorentzian. The deconvolution is easy to solve with basis pursuit: $$ argmin |u|_1 \; subject \; to \; \|l*u - \hat{u_o} \|^2 \leq \mu $$ This ignores $\eta$ as well as our statistical knowledge of $\eta$. Are there ideas on how I can incorporate $\eta$ into my denoising model? Is there a different model I should look into?

edit: looks like I need to set up a MAP estimate for $f = u*\eta + \sigma$. I think I can sort it out when it's just $\eta$ or just $\sigma$.

Answer by rcompton for Best algorithm/software for solving a planar transportation problem ?

2012-03-30T00:16:44Z

How much have you looked into the theory of optimal transport? It's very popular for image warping/registration.

There's codes available to compute the $l1$-optimal transport distance (also referred to as "Earth mover's distance") here: http://ai.stanford.edu/~rubner/emd/default.htm and here: http://www.cs.huji.ac.il/~ofirpele/FastEMD/code/

The $l2$ optimal mass transport problem is quite difficult but can be solved: http://www.springerlink.com/index/40PGJBKDC9V0UH94.pdf

Once it's possible to compute the cost to get between two distributions of points I guess you'll have to optimize to see which distribution is closest to the one you have. Maybe something like: $$ \min_\rho d(\rho_0, \rho)\;subject\;to\;\rho \leq c, \rho \geq 0,\int \rho = 1 $$ where $\rho$ is a probability distribution describing the density of points and $c$ is your threshold. Maybe you can do this with a Lagrange multiplier and gradient descent?

Efficiently computing a few localized eigenvectors

2012-02-11T02:18:54Z

Let $H = \triangle + V(x) : \mathbb{R}^2 \rightarrow \mathbb{R}^2$. I am interested in domain decomposition for an eigenproblem involving $H$.

The lowest 1000 eigenfunctions of $H$, $ \psi_i $, can be partitioned using a region, $\Omega \subset \mathbb{R}^2$, such that each $\psi_i$ localizes either inside of $\Omega$ or outside of $\Omega$. $\Omega$ is not a subspace of $\mathbb{R}^2$ as it may be an oddly shaped region.

Label the inner eigenfunctions $\psi_i^{in}$ and the outer ones $\psi_i^{out}$. There's only about 10 $\psi_i^{in}$s. Given $\Omega$, my goal is to efficiently compute the $\psi_i^{in}$.

One way to find the $\psi_i^{in}$ would be to discretize, compute all 1000 $\psi_i$s, and then partition. This is what I do now (5-point stencil for $\triangle$ on a $10^3 \times 10^3$ grid). The problem is that this requires diagonalizing over a 1000 dimensional space in order to get 10 eigenvectors. It seems like there would be a cheaper way to compute the $\psi_i^{in}$.

Edit: I reposted to http://scicomp.stackexchange.com/questions/1396/efficiently-computing-a-few-localized-eigenvectors#comment2200_1396 and hopefully clarified the problem statement.

Edit I think I can solve this if I can at least figure a way to solve \begin{equation} \max \psi^T H \psi \text{ subject to } P\psi = \psi \text{ and } \psi^T \psi = 1 \end{equation} where $P$ is projection onto the space of functions localized over $\Omega$. My guess is that this will end up looking like power iterations with a projection step built in between matrix applies. If this is doable then something like inverse iteration should be doable which will give me what I want.

Answer by rcompton for How to do (m)Gram-Schmidt orthogonalization with integers ? (real life problem) ("mathematicalized reformulation")

2011-11-07T23:31:58Z

Well, if your microprocessors can handle fixed point arithmetic then here is a matlab commercial that should do it: http://www.mathworks.com/products/fixed/demos.html?file=/products/demos/shipping/fixedpoint/cordicqr_demo.html

Gram-Schmidt is not numerically stable even when you can use floating point so my guess is that you will have many problems if you stay that course.

Definition of spectral gradient

2011-09-20T07:14:56Z

Consider this differential operator $$ \mathcal{H}(\phi(\mathbf{x})) = -\triangle + V(\mathbf{x})H_\epsilon (\phi(\mathbf{x})) $$ where $\mathbf{x} \in \mathbb{R}^2$, $\phi : \mathbb{R}^2 \rightarrow \mathbb{R}^2$, and $H_\epsilon$ is a relaxed version of the Heaviside step function.

I am interested in studying spectra of $\mathcal{H}$ as $\phi$ varies but am not sure about the definitions. Suppose $(v,\lambda)$ is an eigenpair of $\mathcal{H}$, I want to know $\nabla v$ and $\nabla \lambda$.

For example, I compute $\nabla \lambda$ as $$ \mathcal{H}v = \lambda v $$ $$ \implies \nabla (\mathcal{H}v) = \nabla (\lambda v) $$ $$ \implies \nabla \mathcal{H} v + \mathcal{H} \nabla v = \nabla \lambda v + \lambda \nabla v $$ $$ \implies \nabla \lambda = v^T \nabla \mathcal{H} v $$ $$ \implies \nabla \lambda = v^T V(\mathbf{x}) \delta_\epsilon(\phi(\mathbf{x})) v $$

which is weird since $\nabla \lambda$ is a number and I don't know that I computed $\nabla \mathcal{H}$ properly. It's the same way Terry does the computation in the single variable case (cf http://terrytao.wordpress.com/2008/10/28/when-are-eigenvalues-stable/ ).

Any advice? I hope to use $\phi$ in a level set method for domain decomposition.

Edit: Restrict everything to the discrete case (ie $\mathcal{H}$ is a matrix, like in the link above) and I still don't know how to do this. I guess the concise way to state my question is "Given a matrix $\mathcal{H}(\phi)$ does anyone know a formula for the first variation of the eigenvalues or eigenvectors of $\mathcal{H}$?"

Answer by rcompton for Eigenfunctions of random graphs

2011-09-20T07:38:57Z

Laplacian eigenvectors tend to localize over nodes with similar degrees. Some illustrations and exposition begin on page two in this paper: doi:10.1038/nphys1651 with a more detailed explanation here: 10.1103/PhysRevE.77.031102

Operator compression preserving lowest energy eigenspace.

2011-01-28T06:02:57Z

I have a large ($10^6$ by $10^6$) sparse ($0.4$% nonzero) hermitian matrix $H$ arising from the discretization of an elliptic PDE. I would like to approximate $H$ with a smaller matrix $H'$ in such a way that $H'$ and $H$ have nearly identical eigenvectors and eigenvalues for the lowest 5 or 6 eigenvalues.

This could be done if I knew the lowest eigenvectors of $H$ as I could simply restrict $H$ to the space spanned by these values. However I would like to be able to find an approximation before solving eigenvalue problem (in fact my eventual goal is to make solving the eigenvalue problem more efficient).

I know of the recent work on approximating SDD systems with graph sparsification as well as multilevel operator compression. What else is out there?

The application is a full-CI treatment of a multiparticle quantum system.

Answer by rcompton for Identifying poisoned wines

2011-04-02T02:49:24Z

There's a very similar problem in compressed sensing genetic screening for rare alleles (cf http://nar.oxfordjournals.org/content/38/19/e179.full ). The technique almost works here provided we can determine how much poison a rat gets. Seems reasonable, a rat that gets more poison dies faster.

In our problem the idea would be to create a sample for each rat to drink by randomly pooling together wine from many bottles. Specifically, for rat $i$ we draw $A_{ij}$ liters of wine from each of the $j=1,...,N$ bottles where $A_{ij}$ is $\mathcal{N}(0,1)$ distributed. Let $b_i$ denote the amount of poison measured in rat $i$ and let $x_j$ denote the amount of poison in bottle $j$.

This yields the highly underdetermined linear system $$ A\vec{x} = \vec{b} $$ where we know a priori that $\vec{x}$ is sparse. The sparsest solution to this linear system may obtained in polynomial time by solving the convex optimization $$ \min |\vec{x}|_1 \text{ s.t. } A\vec{x}=\vec{b} $$

The number of rats required here is $\mathcal{O}(s \log(N))$ where $s$ is the number of poison bottles.

Answer by rcompton for Remap FFT frequency bin distribution

2011-03-11T00:31:21Z

Pick a large number of points to discretize the frequency domain with. When you have a time signal with less points zero pad until you hit that number. This is sometimes called "spectral interpolation" https://ccrma.stanford.edu/~jos/st/Zero_Padding_Theorem_Spectral.html and does a nice job of interpreting the frequency domain.

Answer by rcompton for Laplace equation with mixed boundary conditions

2010-12-06T12:06:44Z

The solution to the Dirichlet problem is unique by the maximum principle. For the same reason, the solution to the Neumann problem is unique up to a constant.

At worst your problem would fail uniqueness by a constant, but those Dirichlet conditions you've got will prevent this.

This carries over to the discrete case. But note that the discrete case is an easier problem since you can just write out the matrix you need to invert and check if it's invertible.

Answer by rcompton for The Convergence of Jacobi and Gauss-Seidel Iteration

2010-12-03T02:51:34Z

You need to be careful how you define rate of convergence. For Gauss-Seidel and Jacobi you split $A$ and rearrange \begin{eqnarray} Ax &=& b\ M-K &=& b \ \implies x &=& M^{-1}Kx + M^{-1}b\ &\def&Rx + c \end{eqnarray} Giving the iteration $x_{m+1} = Rx_m + c$. We (Demmel's book) define the rate of convergence as the increase in the number of correct decimal places per iteration $$ r = -\log_{10}( \rho(R)) $$ where $\rho(R)$ is the spectral radius of $R$. It can be shown that for $A$ strictly row diagonally dominant that $$ \rho(R_{\text{Gauss}}) \leq \rho(R_{\text{Jacobi}}) < 1 $$ indicating that the rate of convergence for Gauss Seidel is greater than that of Jacobi.

However I have never seen a significant difference in speed between the two methods.

Answer by rcompton for Suggestions for sonifying math

2010-12-02T01:39:38Z

I messed around with this when I was a C++ TA. I came up with "musical sorting algorithms", "musical Gauss Seidel" and a terrible sounding FFT. Details: http://www.math.ucla.edu/~rcompton/art.html

Schrodinger's equation over a randomized grid

2010-11-06T23:26:33Z

I am interested in solutions to $$ \frac{d}{dt} \Psi = -iH \Psi $$ for $H$ hermitian and time independent. This boils down to evaluating $$ \Psi(t) = e^{-iHt}\Psi_0 $$ at points of interest $t_n$. I want to quickly compute $\Psi(t)$ up to a large final time $T$ by spacing the $t_n$ randomly and advancing the propagator over the nonuniform grid.

To be specific, fix $\Delta t$ and draw $t_n$ uniformly at random from $\Delta t, 2\Delta t, ..., N\Delta t$. The gaps $t_{n+1}-t_{n}$ will be distributed geometrically (ie we will have large gaps).

Any ideas how I could compute $e^{-iHt}$ accurately and efficiently when $t$ may be large?

Fast approximations that I know of (eg Strang splitting, Trotter product) are very efficient but I have found them to be horribly inaccurate and unstable for large gaps.

Machine precision accuracy is achieved by expanding $e^{-iHt}$ in Chebyshev polynomials, but the computational cost (measured by counting applications of $H$) scales linearly with the timestep size and we get to $T$ no faster.

I am beginning to think that spectrally accurate approximations of $e^{-iHt}$ provably require linearly more applications of $H$ as $T$ increases. Does anyone know if this is proven?

Thanks for solving my thesis problem,

Ryan

Answer by rcompton for Multivariate Bisection

2010-10-01T02:37:38Z

You might want to consider the vector field

$ \vec{F}(x,y) = (f(x,y), g(x,y)) $

and look for sources and sinks of $\vec{F}$. I think this could be done by recursively dividing up the plane into squares and calculating the winding number of each square. If it is nonzero then you have a critical point within that square (cf Thm 2 this paper) and should divide further.

Answer by rcompton for Is any bias introduced from initial clustering

2010-08-28T00:57:09Z

I've looked at this for PCA-based clustering.

Here is some PC clustering of a simulated dataset. We first draw the clusters as red xs on PCs computed from the entire dataset. A few random points (these are marked blue) are removed from the dataset and the PCs are computed on reduced data. The blue data are then plotted against the reduced PCs and colored green.

I don't have a rigorous bound (I would really like some references) but from what I can tell the most drastic changes are in directions that involve only noise.

Here is some matlab code:

clear all;close all;clc;

%create garbage SNP-style data
m = 500;
n = 25000;
M = zeros(m,n);

%create 2 populations
%one has 0s in the interesting columns
%the other has 2s
interesting_variables = randsample(1:n, .1*n);
M(1:2:m, interesting_variables) = 2.0;

%add sampling noise
sigma = 4.0;
M = M + double(sigma*randn(m,n) > ones(m,n));

%center data
for j=1:n
    M(:,j) = M(:,j) - mean(M(:,j));
end

% plot PCA based clusters, Tygert's code is faster.
tic
[~,~,V] = svds(M,2);
toc


coordsx = zeros(1,m);
coordsy = zeros(1,m);
for i=1:m
    coordsx(i) = M(i,:)*V(:,1);
    coordsy(i) = M(i,:)*V(:,2);
end

figure()
title('strcuture')
plot(coordsx, coordsy, 'rx')


%% 
% Now pick some guys from the dataset to monitor
% identify them as blue triangles in the plot
%

some = m/5;
guys = randsample(1:m, some);

hold on;
coordsx_blue = zeros(1,some);
coordsy_blue = zeros(1,some);
for i=1:some
     coordsx_blue(i) = M(guys(i),:)*V(:,1);
     coordsy_blue(i) = M(guys(i),:)*V(:,2);
end
plot(coordsx_blue, coordsy_blue, 'b<');


%%
% Finally, run PCA again but without the blue guys.
% Observe how this changes the clustering when we plot the blue
% against the PCs found from the reduced data
% green is the new blue.
%

M_small = M;
M_small(guys, :) = [];
[ms ns] = size(M_small);

tic
[~,~,V_small] = svds(M_small,2);
toc

coordsx_newblue = zeros(1,some);
coordsy_newblue = zeros(1,some);
for i=1:some
    coordsx_newblue(i) = M(guys(i),:)*V_small(:,1);
    coordsy_newblue(i) = M(guys(i),:)*V_small(:,2);
end
plot(coordsx_newblue, coordsy_newblue, 'gh');


% draw some lines to see how things moved
for i=1:some
    line([coordsx_blue(i) coordsx_newblue(i)], [coordsy_blue(i) coordsy_newblue(i)],'Color','c');
end

Answer by rcompton for Fast Fourier Transform for Graph Laplacian?

2010-08-12T20:12:23Z

The trick to making the FFT work is factoring out a complex exponential from the sum over odd terms. For this to happen your function needs to be sampled across a uniform grid. Greengard refers to this property as "brittle" (cf math.nyu.edu/faculty/greengar/shortcourse_fmm.pdf ).

When your function is sampled over a nonuniform grid fast multipole methods or Barnes-Hut style algorithms can help.

Answer by rcompton for Switching Research Fields

2010-08-04T03:13:00Z

I like your style. Right now I'm distracting myself from solving a PDE arising in quantum mechanics using techniques from compressed sensing.

Of the 3, compressed sensing is the easiest to pick up as it's more of a new signal processing trick than a fully developed research area. Checking out Candes' ICM talk (http://dsp.rice.edu/cs) as well as Igor's obsessive blog (http://nuit-blanche.blogspot.com/) does not require that much of a time commitment and can be a lot of fun. Finding research problems is another story.

Answer by rcompton for minimize the sum of absolute eigenvalues

2010-07-14T06:25:39Z

Interesting question. Nuclear norm minimization is getting much attention right now as it relates directly to compressed sensing.

Some software for minimization with this constraint that I've used: http://perception.csl.illinois.edu/matrix-rank/sample_code.html

Comment by rcompton

2013-03-06T04:32:12Z

Wait a minute, that paper is about infinite graphs. It makes sense that I can apply it (I have a very large graph) but it's not quite the answer.

Comment by rcompton

2013-02-12T06:03:23Z

I think this will do it. I'll take a closer look at the paper over the next few days, but at a first glance it seems to treat my question exactly.

Comment by rcompton

2013-02-10T02:35:17Z

This is an interesting reference. Are there any English articles related to it?

Comment by rcompton

2013-02-10T02:32:58Z

Oh, right, whoops, this happened on stackexchange to. My graph is very large and has a degree distribution which is (probably) power law. Edited question.

Comment by rcompton

2012-08-23T22:20:51Z

Thanks. I will look into that now.

Comment by rcompton

2012-06-05T23:12:19Z

Also, shouldn't that be argmax?

Comment by rcompton

2012-06-05T23:09:33Z

This is what I think I'll have to do (but I'm not good with Bayesian approaches). It looks similar to Aubert and Aujol's denoising method: www.math.u-bordeaux1.fr/~jaujol/PAPERS/uv.pdf . The sticking point is that I need to minimize something non-convex in $x$. Also, I think we may need to split the $P$ across discretization points and take logs. Otherwise, well need the distribution of $\eta$ over all time which I don't think I can get.

Comment by rcompton

2012-06-04T06:08:28Z

Other researchers have established, using a physical model, that this isn't white noise. A method to solve my problem is still interesting even if we have to ignore that. I'm not sure what I can learn about the correlations in time, the way we can get the density right now is to record only the noise and then use that data to generate noise statistics. We do this at t=0. There might be a way to obtain statistics at future times but I haven't figured out how to do it

Comment by rcompton

2012-06-04T01:38:27Z

We can not approximate it with a well-known distribution. It needs to be recorded each time we run an experiment.

Comment by rcompton

2012-06-04T01:37:17Z

At each point in time we draw $\eta(t_i)$ from a known distribution.

Comment by rcompton

2012-06-03T05:32:13Z

Yes, the $l$ and the $\eta$ hurt the sparsity. The $l$ is no problem since I have an explicit form for it. As for $\eta$, it's noise, all I know is the distribution of $\eta$. I suppose the problem I'm trying to solve is: remove additive noise + remove convolutive noise + remove convolution with $l$.

Comment by rcompton

2012-03-31T01:44:00Z

I would just solve $(I-A)x=b$. What is this iteration you're talking about?

Comment by rcompton

2012-03-30T22:36:33Z

Sure. The book by Villani (who got a fields medal last time around) "Topics in Optimal Transportation" has more information in this field than I can give you.

Comment by rcompton

2012-02-20T07:05:50Z

The setup could be described as domain decomposition for the $\triangle$ operator where $\Omega$ corresponds to one of the domains.

Comment by rcompton

2012-02-18T23:24:21Z

@Squark: the integral is a sum over the indices in $\Omega$. I didn't mention it, but it's mp since $H = \triangle + V(x)$.