Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

For a certain problem I am modelling, I have an MCMC sampler at the moment. It draws samples from the ($n-1$)-dimensional simplex (in this case, from a Dirichlet distribution) and evaluates the respective ratios. The sampler itself works fine, the question comes from transforming the samples:

Say I have $\vec{p}$ from $\mathbb{S}^{n-1}$ embedded in $\mathbb{R}^{n}$ sampled from my Dirichlet distribution. What I am interested in is actually the solution $\vec{f}$ of $\vec{p} = Q^T diag(\vec{p}) \;\vec{f}$, where $Q$ is a stochastic transition matrix, i.e. $\displaystyle \sum_{j=1}^n q_{ij} = 1$ (and is set before). Furthermore, Q is a diagonally dominant matrix, where off-diagonals are usually $10^{-5 \cdot k}$, with $k$ mostly $1, 2$ or $3$ (but can occasionally also go up to say 5 or 6), so precision is important. At the moment I am doing (the deadly sin) of precomputing the inverse of $Q$ and then constantly calculating $\vec{f} = (diag(\vec{p}))^{-1} \; (Q^T)^{-1} \; \vec{p}$.

As I am in the process of rewriting my program in C++ (from MATLAB) I have the following question: what is the best and most accurate way for performing this step, keeping the diagonally dominant matrix in mind? Ultimately, accuracy should be of higher priority, with performance second. As far as I know, Gaussian elimination seems to be the most accurate way to approach this problem, but also scales asymptotically as $\mathcal{O}(n^3)$. For the problem: the dimension $n$ will usually be around 10, and should ultimately never exceed 100. The number of samples $N$ for the sampler will usually be $10^7$.

share
1  
You could precompute a factorization of $Q$ instead of the inverse, which should be numerically more stable and still allow you to evaluate $Q^{-T}p$ in $O(n^2)$ time. –  Rahul Aug 3 '12 at 23:51

This site is currently not accepting new answers.

Browse other questions tagged .