Please clarify- by $[0,1]$, do you mean the closed interval between 0 and 1 in the real numbers, or do you mean the discrete set of values 0 and 1, excluding everything in between?

The fundamental point here is that you don't want to actually compute the Moore-Penrose pseudoinverse, since that matrix will almost certainly be fully dense. Rather, you should consider using an iterative method to compute a minimum norm least squares solution. You could do this by solving a system involving $(A+\frac{ee^{T}}{n})$, or you might want to consider a more straight forward approach. For example, LSQR can be used with a damping factor to minimize a weighted sum of $\| Ax - -b\|^2$ and $\| x \|^2$.

Obviously, they're using different definitions. The notation and conventions associated with the Fourier transform differ between different authors, although it's usually easy to figure out the differences and adjust your results accordingly. A look at the documentation for the R and Mathematica functions should help you figure this out.

It appears that you think you've got {convex} - {concave}, but CVX thinks you've got {convex} + {concave}. I'd carefully check what you're doing and if it still seems write to you, then contact the authors of CVX.

Not really. If r_1 and r_2 are far smaller than r_c, then the additional particle shouldn't have any significant effect on the escape time. The problem only becomes interesting if r_1 is big enough. You also haven't told us what happens if particle 1 and particle 2 interact with each other. Does particle 2 just bounce off particle 1?

You haven't said anything about the relative sizes of the particles and the cylindrical container. Clearly, in one extreme case (P1's radius is as large as the radius of the cylinder so it completely blocks access to part of the cylinder) P1 matters. Just as clearly, P1 doesn't matter in the other extreme case (e.g. the cylinder has a diameter measured in light years and P1 and P2 are the size of small molecules.) You need to tell us something more about what you're actually trying to model...

One particular difficulty with your phase I problem is that many problems have no strictly interior feasible solutions (i.e. the optimal $s$ in your phase I problem is $s^{*}=0$.)

If you'd like an explanation of why Strassen-Winograd and other more sophisticated matrix multiplication algorithms don't work well in practice for matrices of size in the thousands to tens of thousands, I'd be happy to explain that- the basic reason is that Strassen's algorithm trades off matrix multiplications for matrix additions, but matrix additions are relatively more expensive than they should be because matrix addition speed is limited by the poor memory bandwidth of contemporary computers.

Like many other similar computations such as solving a linear programming problem, we're talking about the complexity of obtaining an epsilon approximate solution for some accuracy level. Igor's answer includes a link to a cstheory.stackexchange question whose answers include a link to a STOC paper that appears to discuss how to reduce the eigenvalue decomposition to matrix multiplication. I'm certainly not an expert on this stuff- I'm much more interested in what can be done practically.

you'll need to do this for $c=e_{1}$, $c=e_{2}$, ..., so you'll end up with a block diagonal SDP constraint.

Actually, it depends on the relative size of m and n. If n is much bigger than m, then the multiplication will take the most time. I don't know the size of your m and n, but a simple example with m=100 and n=1000 already showed this effect.