A study of singular spectrum analysis with global optimization techniques

Abstract

Singular spectrum analysis has recently become an attractive tool in a broad range of applications. Its main mechanism of alternating between rank reduction and Hankel projection to produce an approximation to a particular component of the original time series, however, deserves further mathematical justification. One paramount question to ask is how good an approximation that such a straightforward apparatus can provide when comparing to the absolute optimal solution. This paper reexamines this issue by exploiting a natural parametrization of a general Hankel matrix via its Vandermonde factorization. Such a formulation makes it possible to recast the notion of singular spectrum analysis as a semi-linear least squares problem over a compact feasible set, whence global optimization techniques can be employed to find the absolute best approximation. This framework might not be immediately suitable for practical application because global optimization is expectedly more expensive, but it does provide a theoretical baseline for comparison. As such, our empirical results indicate that the simpler SSA algorithm usually is amazingly sufficient as a handy tool for constructing exploratory model. The more complicated global methods could be used as an alternative of rigorous affirmative procedure for verifying or assessing the quality of approximation.

Appendix

In this section we explain how the gradient information can be derived for the global optimization code MultiStart.

Let the objective function (26) be rewritten in the form

$$\begin{aligned} f(\lambda _{1},\ldots ,\lambda _{d}) = \frac{1}{2} \langle \varLambda \mathcal {L}(\varLambda ;\mathbf {z}) - \mathbf {z}, \varLambda \mathcal {L}(\varLambda ;\mathbf {z}) - \mathbf {z} \rangle , \end{aligned}$$

with inner product $\langle \mathbf {p},\mathbf {q} \rangle := \sum _{i} p_{i}\overline{q}_{i}$ for complex vectors. Employing polar coordinates $(\rho _{i},\theta _{i})$ when $\lambda _{i}=\rho _{i}e^{\imath \theta _{i}}$ as variables in our global method, we identify $\varLambda =\varLambda (\lambda _{1},\ldots ,\lambda _{d}) = \varLambda (\rho _{1},\ldots ,\rho _{d}, \theta _{1},\ldots ,\theta _{d})$. We now calculate the gradient of $f$ with respect to the real variables $\rho _{i}$ and $\theta _{i}$.

Using the fact that $\varLambda ^{*} \left( \varLambda \mathcal {L}(\varLambda ;\mathbf {z}) - \mathbf {z}\right) =0$, the action of the Fréchet derivative of $f$ with respect to $\varLambda$ at a complex matrix $H$ is given by

$$\begin{aligned} \frac{\partial f}{\partial \varLambda }.H = \mathfrak {R}\left( \langle H, \left( \varLambda \mathcal {L}(\varLambda ;\mathbf {z}) - \mathbf {z}\right) \mathcal {L}(\varLambda ;\mathbf {z})^{*} \rangle \right) \!, \end{aligned}$$

where $\mathfrak {R}$ stands for the real part of a complex-valued quantity and the same notation $\langle \cdot ,\cdot \rangle$ denotes the generalization to the Frobenius inner product of complex matrices. On the other hand, the action of the Fréchet derivative of $\varLambda$ with respect to variables $\varvec{\rho }=[\rho _{1},\ldots ,\rho _{d}]$ and $\varvec{\theta }=[\theta _{1},\ldots ,\theta _{d}]$ at vectors $\mathbf {h}, \mathbf {k} \in \mathbb {R}^{d}$ can be expressed as

$$\begin{aligned} \frac{\partial \varLambda }{\partial \varvec{\rho }}.\mathbf {h}&= \underbrace{\mathrm {diag}\{0,1,2,\ldots ,2n-2\}}_{\Xi }\varLambda \mathrm {diag}\left\{ \varvec{\rho }\right\} ^{-1} \mathrm {diag}\{\mathbf {h}\}, \\ \frac{\partial \varLambda }{\partial \varvec{\theta }}.\mathbf {k}&= \imath \Xi \varLambda \mathrm {diag}\{\mathbf {k}\}, \end{aligned}$$

respectively. By the chain rule, we obtain the actions of the gradient of $f$ as follows:

$$\begin{aligned} \frac{\partial f}{\partial \varvec{\rho }}.\mathbf {h} = \frac{\partial f}{\partial \varLambda }.\left( \frac{\partial \varLambda }{\partial \varvec{\rho }}.\mathbf {h}\right)&= \mathfrak {R}\left( \langle \Xi \varLambda \mathrm {diag}\left\{ \varvec{\rho }\right\} ^{-1} \mathrm {diag}\{\mathbf {h}\}, \left( \varLambda \mathcal {L}(\varLambda ;\mathbf {z}) - \mathbf {z}\right) \mathcal {L}(\varLambda ;\mathbf {z})^{*} \rangle \right) \!, \\&= \mathfrak {R}\left( \langle \mathrm {diag}\{\mathbf {h}\}, \mathrm {diag}\left\{ \varvec{\rho }\right\} ^{-1}\varLambda ^{*}\Xi \left( \varLambda \mathcal {L}(\varLambda ;\mathbf {z}) - \mathbf {z}\right) \mathcal {L}(\varLambda ;\mathbf {z})^{*} \rangle \right) \!, \\ \frac{\partial f}{\partial \varvec{\theta }}.\mathbf {k} = \frac{\partial f}{\partial \varLambda }.\left( \frac{\partial \varLambda }{\partial \varvec{\theta }}.\mathbf {k}\right)&= \mathfrak {R}\left( \langle \imath \Xi \varLambda \mathrm {diag}\{\mathbf {k}\}, \left( \varLambda \mathcal {L}(\varLambda ;\mathbf {z}) - \mathbf {z}\right) \mathcal {L}(\varLambda ;\mathbf {z})^{*} \rangle \right) \\&= \mathfrak {R}\left( \langle \mathrm {diag}\{\mathbf {k}\}, -\imath \varLambda ^{*}\Xi \left( \varLambda \mathcal {L}(\varLambda ;\mathbf {z}) - \mathbf {z}\right) \mathcal {L}(\varLambda ;\mathbf {z})^{*} \rangle \right) \end{aligned}$$

By the Riesz representation theorem, we see that the gradient of $f$ can be expressed as

$$\begin{aligned} \frac{\partial f}{\partial \varvec{\rho }}&= \mathfrak {R}\left( \mathrm {diag}\left\{ \mathrm {diag}\left\{ \varvec{\rho }\right\} ^{-1}\varLambda ^{*}\Xi \left( \varLambda \mathcal {L}(\varLambda ;\mathbf {z}) - \mathbf {z}\right) \mathcal {L}(\varLambda ;\mathbf {z})^{*} \right\} \right) \!, \\ \frac{\partial f}{\partial \varvec{\theta }}&= \mathfrak {R}\left( \mathrm {diag}\left\{ -\imath \varLambda ^{*}\Xi \left( \varLambda \mathcal {L}(\varLambda ;\mathbf {z}) - \mathbf {z}\right) \mathcal {L}(\varLambda ;\mathbf {z})^{*} \right\} \right) \!. \end{aligned}$$