Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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|cite|improve this question
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

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.