Question in short
Is there an analytical solution for the (non-zero) eigenvalues/pseudo-determinant, i.e., the product of all non-zero eigenvalues, of the real, symmetric, positive semi-definite matrix
$$\mathbf{P} = \mathbf{D}_{m}^{T}\mathbf{D}_{m}$$
where $\mathbf{D}_{m}$ is a finite forward difference matrix?
With $\mathbf{D}_{m}$, the $m$-th order forward difference of a signal $\mathbf{y}$ of size $n$ can be computed as $\mathbf{D}_{m}\mathbf{y}$.
$\mathbf{D}_{m}$ is of shape $\left(n-m\right) \times m$ and $\mathbf{P}$ has $n-m$ non-zero eigenvalues as obvious from the description of $\mathbf{D}_{m}$ below.
Background
I'm trying to implement the Whittaker-Henderson smoother that selects its own smoothing parameter by maximising the marginal likelihood as described by
Biessy G., et al.: Revisiting Whittaker-Henderson Smoothing, 2023, link to arXiv, Chapter (Sections 5 and 5.1, pp. 17-20)
In this specific problem, the pseudo-determinant needs to be computed for the penalty matrix $\mathbf{P}$ as $\ln\left|\mathbf{P}\right|_{+}$. Of course, from a practical point of view, the logarithms of all non-zero eigenvalues should be summed for this. $\mathbf{P}$ is the Gram product of an $m$-th order finite forward difference matrix $\mathbf{D}_{m}$ (the accuracy as named in the Wikipedia articile is 1). So,
$$\mathbf{P} = \mathbf{D}_{m}^{T}\mathbf{D}_{m}$$
where, e.g.,
$$\mathbf{D}_{1} = \begin{pmatrix} -1&1&0&0&0&0\\ 0&-1&1&0&0&0\\ 0&0&-1&1&0&0\\ 0&0&0&-1&1&0\\ 0&0&0&0&-1&1 \end{pmatrix}$$
$$\mathbf{D}_{2} = \begin{pmatrix} 1&-2&1&0&0&0\\ 0&1&-2&1&0&0\\ 0&0&1&-2&1&0\\ 0&0&0&1&-2&1\\ \end{pmatrix}$$
for obtaining the finite forward differences of a signal $\mathbf{y}$ of size $n=6$. Only rows with full coefficients are used, leaving a shape of $\left(n-m\right) \times n$.
In this context, $m$ can go much higher, e.g., for spatially adaptive smoothing that involves the 6-th order difference for obtaining a smooth estimate of the 4-th order derivative via the Whittaker-Henderson smoother.
Given all this, it is easy to obtain that $\mathbf{P}$ is a real, symmetric, positive semi-definite matrix with $n-m$ non-zero eigenvalues.
Problem
I'm struggling to find an efficient and accurate way to obtain the log-pseudo-determinant $\ln\left|\mathbf{P}\right|_{+}$, especially for higher orders of $m$ and large signal sizes $n$ because the computation time increases dramatically for $n > 1000$ and for $m > 4$ the ratio between the smallest non-zero eigenvalue $\lambda_{\min.}^{+}$ and the largest eigenvalue $\lambda_{\max.}$ exceeds $10^{15}$ already for $n < 1000$, thereby causing $\lambda_{\min.}^{+}$ and its closest neighbours to be possibly slightly negative since they are computed from 64-bit-float numbers. This prohibits the computation of $\ln\left|\mathbf{P}\right|_{+}$ since then it is impossible to sum the logarithms of the eigenvalues if only one goes negative as its presence in $\ln\left|\mathbf{P}\right|_{+}$ is crucial despite its magnitude.
I was thus thinking that an analytical solution (or even a crude approximation only) for the eigenvalues/log-pseudo-determinant would
- be cheaper to compute
- not suffer as much from floating point inaccuracies since the eigenvalues are not originating form a cascade of operations on the same matrix.
What I tried so far
First, I was tackling the problem numerically using a dense eigenvalue solver which seems to suffer less from floating point precision, but cannot be scaled to large $n > 2500$ as it takes more than 10 seconds to compute (Python NumPy Hermitian Eigensolver which relies on LAPACK).
Afterwards, I tried another dense solver which is faster, but apparently at the expense of precision causing negative eigenvalues for $m=4$ at $n > 1000$ (Python SciPy Hermitian Eigensolver with also relies on LAPACK, called with “evr” driver).
From looking at the R-implementation of the original publication named above, it is apparent that a dense solver was used without any special measures against floating point imprecision.
Then - since $\mathbf{P}$ is banded with $m$ off-diagonals on each side - I tried to use an Eigensolver for Hermitian Banded Matrices which relies on a divide and conquer algorithm. This proved to be more than twice as fast as the dense solvers, but also gave negative eigenvalues for $m=4$ and $n > 800$ (Python SciPy Hermitian Banded Eigensolver).
Finally, I tried to follow the Wikipedia article on the eigendecomposition of the second order derivative with different boundary conditions because I found that this is closest to the problem at hand from all online sources I dug into. That's where I'm stuck because I don't find a way to re-phrase $\mathbf{P}$ with arbitrary $m$ to match this approach - if this is possible at all.
Alternative Solutions
As mentioned above, a crude approximation of the eigenvalues or the log-pseudo-determinant would probably be a huge step forward.
Otherwise, it could maybe also be beneficial to re-think the boundary conditions of $\mathbf{D}_{m}$ as long as they are not periodic. I was thinking about a mirrored ($d c b \left| a b c d e f g \right| f e d$) or repeating boundary ($a a a \left| a b c d e f g \right| g g g$).
Related Questions
There is this question that deals with higher order finite difference matrix eigendecomposition, but this refers to different accuracies and not orders if the naming from the Wikipedia article on finite differences is applied.
Thanks very much for any answers and comments!
Note: I'm familiar with basic linear algebra, but not a professional and thus not that used to all mathematical abbreviations and symbols