This is a follow-up to estimate precision matrix with given spatial sparsity pattern, expanding on the second part of that question and formulating more precisely using material from the answer by jlewk.

I have a matrix of measurements $\Xi$, where each of the $n$ rows comprises a measurement of $p$ variables. I assume that the variables are normally distributed with expectation zero and that the measurements are i.i.d., but the variables can be correlated and of different variance. It holds $p \gg n$, and $p$ can be on the order of several 100,000.

I would like to estimate a $p \times p$ whitening transformation $W$ such that the whitened data are distributed as $$ \Xi W \sim \mathcal{MN} (0, I_n, I_p), $$ where $\mathcal{MN}$ denotes the matrix normal distribution. Maximum likelihood estimation leads to the function $$ \DeclareMathOperator{\tr}{tr} \newcommand{\dt}[1]{\left\lvert#1\right\rvert} \tr(\Xi' \Xi W W') - n \log \dt{W W'} $$ to be minimized. Without further constraints, this leads to $(W W')^{-1} = \frac1n \Xi' \Xi$, but because $p > n$, $\Xi' \Xi$ is not invertible. Note that $P = W W'$ is the inverse covariance or precision.

To circumvent this estimation problem, I want to impose a sparsity constraint on $W$.

The different variables $\xi_i$ (columns of $\Xi$) are measured at different positions in space $(x_i, y_i, z_i)$, and I know that correlations between them are mainly due a local signal mixture, which would make the precision matrix sparse. Moreover, for technical reasons I would prefer each whitened variable (column of $\Xi W$) to derive only from a small set of original variables near to each other. This motivates the sparsity constraint $$ W_{ij} = 0 \quad \text{for} \quad d_{ij} > r, $$ where $$ d_{ij} = \sqrt{ (x_i - x_j)^2 + (y_i - y_j)^2 + (z_i - z_j)^2 } $$ is the distance between variable measurement positions and $r$ is an upper limit on this distance.

I'm aware that the index $j$ in $W_{ij}$ refers to a whitened variable which strictly speaking does not have a position, but the sparsity constraint with respect to $W$ still makes sense because the implied precision $W W'$ fulfills approximately the same sparsity constraint with a radius of $2 r$.

My question is: How can I efficiently calculate the $W$ that minimizes the function under the constraint?

I can of course use a generic optimization algorithm and enforce the constraint by using parameters only for the non-zero elements of $W$ – but I would hope that there is a better way.

jlewk answered the same question with respect to the precision $P$ by pointing to variants of the graphical lasso where a penalty is applied not uniformly to all elements but element-specific, $$ \sum_{i \neq j} \Lambda_{ij} \dt{P_{ij}}, $$ so that a given sparsity pattern can be approximately enforced by choosing $\Lambda$ appropriately; and to implementations of these variants in different software packages. However, in this new question penalties would have to be applied to the elements of $W$, and the minimization is complicated by $P$ being quadratic in $W$.

In detail:

– Has this problem, "sparse localized whitening estimation", been considered in the literature before?

– Is the mathematical form of the constrained optimization problem outlined above known, and are there standard approaches how to solve it?

– Are there specific implementations, or known implementation strategies, that make this estimation process particularly efficient in light of my large $p$?


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.