Loading repository data…
Loading repository data…
sliaghat / repository
This repository contains a CUDA/C++ implementation of a distinct algorithmic approach for finding the linear recurrence relation (with least possible order) of a real-valued, noisy 1D sequence in near-linear time.
A transparent discovery signal based on current public GitHub metadata.
This score does not audit code, security, maintainers, documentation quality, or suitability. Verify the repository and its current documentation before adoption.
This repository contains a CUDA/C++ implementation of a distinct new algorithmic approach (to the best of my knowledge) for finding the linear recurrence relation (with least possible order) of a real-valued, noisy 1D sequence that may satisfy it. This project features a highly optimized, GPU-accelerated solver that exploits the properties of this problem to reduce the overall time complexity to $O(I \cdot L \log^2 L)$ (where $L$ is the dimension of the Hankel matrix constructed from the given sequence and $I$ is the number of LSQR iterations for solving linear least square as a subproblem).
Given a noisy sequence $S = (s_0, s_1, \dots, s_{2L-1})$, the goal is to find the minimum integer $k$ (also called the linear complexity) and the corresponding coefficients $x = (x_0, x_1, \dots, x_{k-1})$ such that any sequence element $s_{i}$ can be approximated by a linear combination of its $k$ previous terms:
$$ s_{k+i} \approx \sum_{j=0}^{k-1} s_{i+j} x_j $$
The given sequence can conceptually be mapped to a square Hankel matrix $H$ of size $L \times L$:
$$ H = \begin{bmatrix} s_0 & s_1 & s_2 & \dots & s_{L-1} \ s_1 & s_2 & s_3 & \dots & s_L \ \vdots & \vdots & \vdots & \ddots & \vdots \ s_{L-1} & s_L & s_{L+1} & \dots & s_{2L-2} \end{bmatrix} $$
Finding the linear complexity $k$ of the sequence is equivalent to finding the rank of $H$. In other words, we have to find the first row in $H$ that is linearly dependent on its preceding rows.
In order to check if the $k$-th row is dependent on rows $0$ through $k-1$, we can formulate a linear least-squares problem as follows:
$$ A x \approx b $$
Where:
$$ A = \begin{bmatrix} s_0 & s_1 & \dots & s_{k-1} \ s_1 & s_2 & \dots & s_k \ \vdots & \vdots & \ddots & \vdots \ s_{L-1} & s_L & \dots & s_{k+L-2} \end{bmatrix} $$
We then solve $\min_x ||Ax - b||_2$. If the relative residual $\frac{||Ax - b||_2}{||b||_2}$ is below a defined threshold $\epsilon$, the row can be considered linearly dependent.
In a Hankel matrix which represents a linear recurrence relation, if row $k$ is linearly dependent on its preceding rows, then row $k+1$ is guaranteed to be dependent on its preceding rows as well. Because of this property that guarantees a single transition point from independent to dependent, we only need to find the first linearly dependent row, which can be done via binary search. Here a binary search can correctly find the transition point (the linear complexity $k$) in $O(\log L)$ checks, bypassing the need for a linear $O(L)$ search over the rows.
The LSQR algorithm only requires computing two operations at each of its iterations: $y \leftarrow y + Ax$ and $x \leftarrow x + A^Ty$. Given the structure of Hankel matrices, the matrix-vector multiplication $Ax$ is exactly the discrete cross-correlation between the sequence $S$ and the vector $x$:
$$ (A x)i = \sum{j=0}^{k-1} s_{i+j} x_j $$
We know that by the Convolution Theorem, cross-correlation can be computed in the frequency domain using the complex conjugate of the Fourier-transformed vector:
$$ A x = \mathcal{F}^{-1} \Big( \mathcal{F}(S) \odot \overline{\mathcal{F}(x)} \Big) $$
Where $\odot$ is element-wise complex multiplication and $\overline{Z}$ is the complex conjugate. This exact same logic can be applied to $A^T y$:
$$ (A^T y)j = \sum{i=0}^{L-1} s_{i+j} y_i \implies A^T y = \mathcal{F}^{-1} \Big( \mathcal{F}(S) \odot \overline{\mathcal{F}(y)} \Big) $$
Using FFT, this reduces the time complexity of each LSQR iteration from $O(L \cdot k)$ to almost $O(L \log L)$.
Because $H$ is only defined by $S$, we can bypass building $H$. So space complexity drops: