Higher-order character decomposition: char module (Checked: D)#

NOTE: This tutorial assumes familiarity with the concepts explained in Conceptual and mathematical background of higher-order Fourier analysis and that the HoFa package has already been installed, see Installation.

In this notebook we are going to learn how the HoFa package isolates the higher-order components of a function.

Important considerations#

As we recall from Conceptual and mathematical background of higher-order Fourier analysis, the question of what the higher-order components of a function are is subtle, because in higher-order Fourier analysis (when we go beyond classical Fourier analysis) we lose the concept of exactness and replace it with the concept of sparsity and efficiency.

This means that, a priori, we cannot expect an exact decomposition of a function \(f\) into \(f=f_e+\sum g_i\) where \(f_e\) is an error term and the \(g_i\) are the unique higher-order characters. Nevertheless, there exist approximate concepts that are reasonable.

REMARK: Remember that Discrete Fourier Transform gives a unique decomposition for functions on (say) \(\mathbb{Z}/n\mathbb{Z}\) and thus, any higher-order Fourier character can be decomposed uniquely as some linear combination of Fourier characters. However, in general this decomposition may not be sparse or efficient. That is why we need higher-order characters.

A natural generalization#

Throughout this package, we have talked about a natural generalization of simple Fourier characters such as

\[ x\mapsto \exp(2 \pi i x \xi)\]

are quadratic Fourier characters, such as

\[ x\mapsto \exp(2 \pi i (x^2 \xi+\eta x)).\]

It turns out that this intuition is correct and such functions can be regarded as higher-order Fourier characters (or \(k\)-th order Fourier characters if we want to emphasize the particular order \(k\)). In the case of the \(U^3\)-norm, we will denote them by quadratic Fourier characters.

Let us see how the HoFa package can be used to find such a decomposition.

A simple example#

Let us begin by importing the necessary libraries. Recall that we are assuming that the HoFa package is already installed, see Installation. For a quick installation, run the following command:

(hofavenv) $ pip install hofa

import numpy as np
import matplotlib.pyplot as plt
np.set_printoptions(precision=5)

import hofa.char as char
import hofa.rgz as rgz
import hofa.cdf as cdf

Consider a single quadratic Fourier character \(g(x):=e^{2\pi i x^2/n}\), for \(x=0,\ldots,n-1\). As usual, imagine that we perturb it with some random noise \(f_r(x)\sim \text{Unif}[-1,1)\) in a way that we can only see the convex combination \(f:=tf_r+(1-t)g\) for \(t=0.3\).

n = 101
t = 0.4
np.random.seed(123)

x = np.arange(n)
g = np.exp(2*1j*np.pi*x**2/n)

f = np.random.uniform(-1, 1, size=len(x))*t+(1-t)*g

The hofa.char module contains a method spechoft (Spectral higher-order Fourier transform) which allows us to recover the underlying structured component \((1-t)g\). It returns several data relating to the process of finding those characters, but the most important part is the attribute higher_order_char, which is a numpy.ndarray of shape (f.shape,m) where m represents the number of components (higher-order characters) found.

If we represent the result as a family of vectors \(v_i\), the higher-order character decomposition of \(f\) is as follows:

\[ f = f_e+\sum_i \langle f,v_i\rangle v_i \text{ where } f_e \text{ is the error term.}\]

In our case, since we know that the function \(f\) in our example will contain exactly one higher-order character, we simply define v below as such an element.

v = char.spechoft(f, rng=123).higher_order_char[...,0]

coeff_proj = np.mean(f*v.conjugate())

error = np.mean(np.abs(g*(1-t)-coeff_proj*v))
print(r'L1 error between the known quadratic Fourier character (1-t)g and <f,v>v : '+str(error))

L1 error between the known quadratic Fourier character (1-t)g and <f,v>v : 0.012892689137707978

As we can see, the error is very small and, up to the scaling factor \((1-t)\), we are able to recover the character \(g\).

Let us plot below the (real part of the) original function \((1-t)g\), the (real part of the) noisy version \(f\), and the (real part of the) reconstruction which comes from the projection of \(f\) to the (sole) eigenvector returned by spechoft.

fig, ax = plt.subplots(figsize=(12, 5))

# Plot original (1-t)*f_s in light green
ax.plot(x, np.real(g*(1-t)), linestyle='-', color='lightgreen', marker='o', markersize=6, linewidth=1, alpha=0.6, label=r'$Re((1-t)g)$')

# Plot modified function in red
ax.plot(x, np.real(f), linestyle='-', color='red', marker='o', markersize=6, linewidth=1, alpha=0.6, label=r'$Re(f) = Re(tf_r+(1-t)g)$')

# Plot recovered function in blue
ax.plot(x, np.real(coeff_proj*v), linestyle='-', color='blue', marker='o', markersize=3, linewidth=1, alpha=0.6, label=r'$Re(\langle f,v\rangle v)$ where v is the eigenvector found')

# Customize the plot
ax.set_xlabel("x")
ax.set_title(f"t (Amount of noise): {t:.2f}")
ax.set_ylim(-1,1)

# Legend below the plot
fig.subplots_adjust(bottom=0.25)  # Adjust layout to fit legend below
ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.15), ncol=2, frameon=False)

# Show the plot
plt.show()
../_images/05f988047861b8de75415dbcd8a764198ccc7b45dd1352c438e5383a22d17ab8.png

However, this example may not seem very impressive because, after all, we already know a method from the hofa package that can deal with it, hofa.rgz.regularize.

The real challenge is when there are several higher-order Fourier characters that cannot be distinguished by simply applying hofa.rgz.regularize. In the next example we will see how the package hofa.char contains a method that can.

A more interesting and important example#

Now, let \(g(x):=\sin(2\pi x^2/n)\) and, as above, let \(f_r(x)\sim \text{Unif}[-1,1)\). Again, assume that we are presented with a noisy version \(f:=tf_r+(1-t)g\) for some \(t\).

Returning to our intuition with quadratic Fourier characters, it is reasonable to aim to find that \(g\) can be decomposed as

\[ g(x) = \sin(2\pi x^2/n) = \frac{\exp(2\pi i x^2/n)-\exp(-2\pi i x^2/n)}{2i}. \]

Or, equivalently, that the quadratic Fourier characters that are hidden inside \(f\) are \((1-t)\exp(2\pi i x^2/n)/2i\) and \(-(1-t)\exp(-2\pi i x^2/n)/2i\).

n = 1501
t = 0.3

x = np.arange(n)
g = np.sin(2*np.pi*x**2/n)

# We fix a seed for reproducibility
np.random.seed(28)
f = np.random.uniform(-1, 1, size=len(x))*t+(1-t)*g

The method char.spechoft is precisely designed to find such components

# As usual, we fix an rng for reproducibility
rnd_sep_result = char.spechoft(f, rng=28)

eichars = rnd_sep_result.higher_order_char

In the attribute higher_order_char we have an array containing the quadratic Fourier characters that make up \((1-t)g\), i.e. the structured part of \(f\). In this case (check it in your own notebook!), we find exactly two characters that essentially correspond to the two quadratic Fourier characters that make up \(g\).

Below, we plot the different elements of this simple example. For illustration, we plot first the functions \(f\) and \((1-t)g\), both of which are real-valued, and for the quadratic Fourier characters and the quadratic components found, we plot both their real and imaginary parts. Note that, analogously to what happens in the case of the standard Fourier transform (using complex exponentials as Fourier characters), even if the function that we want to decompose is real-valued, e.g. \(f\) and \((1-t)g\), its decomposition may require complex numbers.

import itertools
import numpy as np
import matplotlib.pyplot as plt

number_of_first_elems_to_plot = 500

v_1 = eichars[:, -1]
coeff_proj_1 = np.mean(f * v_1.conjugate())

v_2 = eichars[:, -2]
coeff_proj_2 = np.mean(f * v_2.conjugate())

# Restrict to plotting window
x_plot = x[:number_of_first_elems_to_plot]

# Target ("green") complex curves
target_1_complex = (1 - t) * np.exp(2 * np.pi * 1j * x_plot**2 / n) / (2*1j)
target_2_complex = -(1 - t) * np.exp(-2 * np.pi * 1j * x_plot**2 / n) / (2*1j)

targets_complex = [target_1_complex, target_2_complex]

# Recovered ("blue") complex curves
blue_1_complex = coeff_proj_1 * v_1[:number_of_first_elems_to_plot]
blue_2_complex = coeff_proj_2 * v_2[:number_of_first_elems_to_plot]

blues_complex = [blue_1_complex, blue_2_complex]

# --- L^2 (squared) error comparison over all possible assignments ---

best_perm = None
best_err = np.inf

for perm in itertools.permutations(range(2)):
    total_err = sum(
        np.sum(np.abs(blues_complex[perm[i]] - targets_complex[i])**2)
        for i in range(2)
    )
    if total_err < best_err:
        best_err = total_err
        best_perm = perm

# Reorder blues according to best assignment
matched_blues_complex = [blues_complex[best_perm[i]] for i in range(2)]

# Labels
target_labels_real = [
    r'$\mathrm{Re}\!\left((1-t)\exp(2\pi i x^2/n)/2i\right)$',
    r'$\mathrm{Re}\!\left(-(1-t)\exp(-2\pi ix^2/n)/2i\right)$',
]

target_labels_imag = [
    r'$\mathrm{Im}\!\left((1-t)\exp(2\pi i x^2/n)/2i\right)$',
    r'$\mathrm{Im}\!\left(-(1-t)\exp(-2\pi ix^2/n)/2i\right)$',
]

blue_labels_real = [
    r"real part of matched quadratic Fourier character 1",
    r"real part of matched quadratic Fourier character 2",
]

blue_labels_imag = [
    r"imaginary part of matched quadratic Fourier character 1",
    r"imaginary part of matched quadratic Fourier character 2",
]

# --- Plotting ---

# 4 rows x 2 columns:
# top row spans both columns,
# then 3 rows of 2 plots each
fig = plt.figure(figsize=(15, 16))
gs = fig.add_gridspec(4, 2, height_ratios=[1, 1, 1, 1], width_ratios=[1, 1])

# Top plot
ax_top = fig.add_subplot(gs[0, :])
ax_top.plot(
    x_plot,
    f[:number_of_first_elems_to_plot],
    linestyle="-",
    color="red",
    marker="o",
    markersize=5,
    linewidth=1,
    alpha=0.6,
    label=r"input $f$ (noisy version given to the algorithm)",
)
ax_top.plot(
    x_plot,
    (1 - t) * g[:number_of_first_elems_to_plot] / 3,
    linestyle="-",
    color="green",
    marker="o",
    markersize=5,
    linewidth=1,
    alpha=0.6,
    label=r"scaled $g$ (quadratic truth)",
)
ax_top.set_xlabel("x")
ax_top.set_title(f"t (Amount of noise): {t:.2f}, n (total length): {n}")
ax_top.set_ylim(-1, 1)
ax_top.legend(loc="upper center", bbox_to_anchor=(0.5, -0.15), ncol=2, frameon=False)

# Bottom 3x2 comparison plots
for i in range(2):
    target = targets_complex[i]
    blue = matched_blues_complex[i]

    # Left: real parts
    ax_real = fig.add_subplot(gs[i + 1, 0])
    ax_real.plot(
        x_plot,
        np.real(target),
        linestyle="-",
        color="lightgreen",
        marker="o",
        markersize=5,
        linewidth=1,
        alpha=0.6,
        label=target_labels_real[i],
    )
    ax_real.plot(
        x_plot,
        np.real(blue),
        linestyle="-",
        color="blue",
        marker="o",
        markersize=3,
        linewidth=1,
        alpha=0.6,
        label=blue_labels_real[i],
    )
    ax_real.set_xlabel("x")
    ax_real.set_ylim(-1, 1)
    ax_real.legend(loc="upper center", bbox_to_anchor=(0.5, -0.15), ncol=1, frameon=False)

    # Right: imaginary parts
    ax_imag = fig.add_subplot(gs[i + 1, 1])
    ax_imag.plot(
        x_plot,
        np.imag(target),
        linestyle="-",
        color="lightgreen",
        marker="o",
        markersize=5,
        linewidth=1,
        alpha=0.6,
        label=target_labels_imag[i],
    )
    ax_imag.plot(
        x_plot,
        np.imag(blue),
        linestyle="-",
        color="blue",
        marker="o",
        markersize=3,
        linewidth=1,
        alpha=0.6,
        label=blue_labels_imag[i],
    )
    ax_imag.set_xlabel("x")
    ax_imag.set_ylim(-1, 1)
    ax_imag.legend(loc="upper center", bbox_to_anchor=(0.5, -0.15), ncol=1, frameon=False)

fig.tight_layout()
plt.show()
../_images/65e190eafd6a530f2e1ddd9d3cf7e8df14c937681d8cfd5597ac1467b632fe18.png

How the algorithm works#

Recall from Tutorial - Denoising (regularization), rgz and cdf that the regularization process implemented in hofa.rgz.regularize involves, for a given function \(f\), finding some significant, i.e. large, eigenvalues (of some matrix \(\mathcal{K}(f\otimes \overline{f})\)) and then projecting \(f\) onto the corresponding eigenvectors. it is then natural to think that those eigenvectors should be the higher-order characters that underlie our function \(f\). However, this is not always the case, as we will see in the following example phrased for clarity in classical Fourier analysis.

Why separation of eigenvalues is needed#

Let \(n\) be a positive integer and let \(f:\{0,\ldots,n\}\to \mathbb{C}\) be the function defined by \(f(x) = \cos(2\pi x/n)\). Suppose we want to find the Fourier components of this function using only our regularization algorithm and linear algebra. That is, we do not have access to the knowledge of what the Fourier characters are on the interval \(\{0,\ldots,n\}\).

The way in which hofa.rgz.regularize would operate (importantly, when used with order=1) is to create a matrix \(\mathcal{K}(f\otimes \overline{f})\) which is equal to

\[ \tfrac{1}{4}\chi_1\otimes \overline{\chi_1} +\tfrac{1}{4}\chi_2\otimes \overline{\chi_2}\]

where \(\chi_1(x) = \exp(2\pi i x/n)\) and \(\chi_2(x) = \exp(-2\pi i x/n)\).

NOTE: The interested reader may try to derive this result. The idea is to apply to each diagonal the average operator. That is, replace each \(Z\)-diagonal by the constant value equal to its average.

The critical point here is that, from the perspective of linear algebra, the matrix \(\mathcal{K}(f\otimes \overline{f})\) is a rank 2 matrix with 2 eigenvectors of equal eigenvalue \(\tfrac{1}{4}\). Hence, if we compute its spectral decomposition, again only looking at it from the perspective of linear algebra, this matrix has one eigenspace of dimension 2 associated with the eigenvalue 1/4 and another eigenspace of dimension n-2 and eigenvalue 0.

But this poses a problem for finding the individual Fourier characters \(\chi_1,\chi_2\), because from the perspective of linear algebra we see only a single eigenspace where, a priori, we cannot separate properly these two characters. In practical terms, if we tried to use a linear algebra module to find the eigenvalues and eigenvectors of that matrix we would get the correct eigenvalues but as eigenvectors we would get vectors of the form

\[ a \chi_1+b\chi_2 \text{ and } b \chi_1-a\chi_2 \text{ where } |a|^2+|b|^2 = 1.\]

A solution: take a random sample and repeat the algorithm#

Starting from the previous example, once we have found the eigenspace generated by \(\chi_1,\chi_2\) we could randomly sample a vector of modulus 1 from such eigenspace. Note that such a random vector will have the form \(c \chi_1+d\chi_2 \) for \( |c|^2+|d|^2 = 1\). The critical observation here is that with high probability (in fact, in this case, with probability 1) we get that \(|c|\not=|d|\). Hence, if we apply again the method hofa.rgz.regularize (with order=1) to the function \(c \chi_1+d\chi_2 \) then the corresponding matrix will be precisely

\[ |c|^2\chi_1\otimes \overline{\chi_1} +|d|^2\chi_2\otimes \overline{\chi_2}.\]

Therefore, we are done! Now we have two eigenspaces with different eigenvalues, and thus any linear algebra method will correctly give as eigenvectors \(\chi_1\) and \(\chi_2\).

How this translates for the higher-order cases#

As we have seen before, exact results in classical Fourier analysis translate to approximate results in higher-order Fourier analysis. Thus, concepts such as having different eigenvalues become having eigenvalues sufficiently separated. In that spirit, let us present a minimalistic pseudocode of how the function hofa.char.spechoft works.

Minimalistic pseudocode of the spectral higher-order Fourier transform

Input: Function to process: \(f: Z \to \mathbb{C}\) and order of regularization \(k\).

  1. Initial regularization:

    1. \(f_{\text{reg,init}} \leftarrow \text{regularize}(f,k)\)

    2. \(\rho_{\text{init}} \leftarrow \text{initial\_threshold}(f, f_{\text{reg,init}})\)

    3. \(N_{\text{target}} \leftarrow \# \{\text{eigenvalues of } f_{\text{reg,init}} \text{ greater than } \rho_{\text{init}} \}\)

  2. Initialization for the loop:

    1. \(h \leftarrow f_{\text{reg,init}}\quad\)

    2. \(\text{are\_separated} \leftarrow \text{False}\)

  3. While \(\text{are\_separated} == \text{False} \):

    1. \(f_{\text{reg,iter}} \leftarrow \text{regularize}(h,k)\)

    2. \(\rho_{\text{iter}} \leftarrow \text{iteration\_threshold}(h, f_{\text{reg,iter}})\)

    3. \(N_{\text{found}} \leftarrow \# \{\text{eigenvalues of } f_{\text{reg,iter}} \text{ greater than } \rho_{\text{iter}} \}\)

    4. If \(N_{\text{found}} == N_{\text{target}}\):

      • \(\text{are\_separated} \leftarrow \text{are\_eigenvalues\_sufficiently\_separated}(\text{eigenvalues of }f_{\text{reg,iter}})\)

    5. If \(\text{are\_separated} == \text{False}\):

      • \(h \leftarrow \text{random\_unit\_vector}(\{\text{eigenvalues of } f_{\text{reg,init}} \text{ greater than } \rho_{\text{init}} \})\)

Output: Eigenvectors of \(f_{\text{reg,iter}}\)

While the idea behind this pseudocode is the core of hofa.char.spechoft, there are many technical improvements to consider. For instance, clearly we do not want a potentially infinite loop to run and the thresholds in the pseudocode have to be computed somehow. Moreover, there is a conceptual improvement that we can do: instead of randomly sampling repeatedly until we find a function such that all eigenvalues are sufficiently separated, we could keep track of separated eigenvalues, subtract from \(f\) the contribution of its corresponding eigenvectors (i.e. for each eigenvector \(v\) whose eigenvalue is separated from the rest we replace \(f\) by \(f-\langle f,v\rangle v\) ), keep track of those, and then repeat the algorithm trying to find eventually less and less components.

Algorithm refinements#

The pseudocode above captures the essential idea behind the spectral higher-order Fourier transform implemented in hofa.char.spechoft. Let’s expand on this core concept and discuss the important refinements that make the algorithm practical and efficient.

Let us discuss the possible refinements of the pseudocode Minimalistic pseudocode of the spectral higher-order Fourier transform:

  • In step 1.A recall that hofa.rgz.regularize offers the possibility of choosing different forms of regularization depending on certain hofa.rgz.LayerRegularizer. Similarly, hofa.char.spechoft offers the possibility of passing a list of hofa.rgz.LayerRegularizer to perform the initial regularization.

  • In step 1.B a threshold must be computed, and hofa.char.spechoft is designed in a way that the user may define how to compute this value.

  • In step 3. the loop may run indefinitely. Hence, hofa.char.spechoft always requires a maximum number of iterations that may depend on \(f\) and on the initial regularization \(f_{\text{reg,init}}\).

  • In step 3.A we need to perform a second regularization, this time of \(h\) (which is first \(f_{\text{reg,init}}\) and later changes). To perform this regularization we may use the same regularizers that we used in step 1.A but we may want to use a different one. The function hofa.char.spechoft allows a different regularization procedure given by the user.

  • In step 3.B we have the same issue as in step 1.B, there is a threshold that must be computed and hofa.char.spechoft allows the user to choose how to do it.

  • In step 3.D we can do two things:

    • Try to do as in the pseudocode and declare success for the hofa.char.spechoft only if we have found some \(h\) such that all its eigenvalues are sufficiently separated.

    • Perform an iterative search, where for any given \(h\), if we let \(v_1,\ldots,v_m\) be the eigenvectors such that their corresponding eigenvalues are sufficiently separated, record those values, and then let \(f\leftarrow f-\sum_{i=1}^m \langle f,v \rangle v\).

  • In step 3.B we have to compute some kind of threshold to declare when some eigenvectors are sufficiently separated.

The Role of Randomization#

Randomization plays a crucial role in the algorithm, as it allows us to find higher-order components (characters) of a function analogously to what happened in the example for classical Fourier analysis.

Convergence and Termination#

The algorithm terminates when one of the following conditions is met:

  • Success: All target higher-order characters have been successfully separated either iteratively or all at once.

  • Iteration Limit: The maximum number of iterations is reached without success.

The iterative approach, where successfully separated components are recorded and removed from further consideration, significantly improves the algorithm’s efficiency.

The Random Separator Policy#

The algorithm’s flexibility is implemented via a policy object that encapsulates the strategy for determining Algorithmic refinements. This is where the hofa.char.RandomSeparator abstract base class comes into play. Its purpose is to determine:

  • Initial Threshold: How to compute the starting minimum threshold for identifying potential higher-order components.

  • Iteration Threshold: How to compute the minimum thresholds during the iterative refinement process.

  • Separation Gap: How to determine when eigenvalues are sufficiently separated.

  • Iteration Control: How many iterations to perform

  • Iterative search: whether to remove already found characters or to try to find a candidate random function such that all its significant eigenvalues are already sufficiently separated.

This separation of concerns allows the core algorithm in spechoft to remain generic while enabling customization through different separator implementations.

The class hofa.char.RandomSeparator#

The RandomSeparator abstract base class defines the interface for the algorithm. This interface is designed to encapsulate the core functionality needed for hofa.char.spechoft, providing a flexible framework that can be adapted to various specific implementations. Let us see the definition in Python of such abstract class.

# In hofa.char
class RandomSeparator(ABC):
    @abstractmethod
    def initial_threshold(self, f: np.ndarray, reg_res: rgz.RegularizationResult) -> float:
        """Compute the initial threshold for eigenvalue separation."""

    @abstractmethod
    def iteration_threshold(self, f: np.ndarray, reg_res: rgz.RegularizationResult,
                        initial_f: np.ndarray, initial_reg_res: rgz.RegularizationResult) -> float:
        """Compute the threshold for each iteration of the separation process."""

    @abstractmethod
    def separation_gap(self, f: np.ndarray, reg_res: rgz.RegularizationResult,
                   initial_f: np.ndarray, initial_reg_res: rgz.RegularizationResult) -> float:
        """Compute the gap required between eigenvalues for successful separation."""

    def get_rng(self) -> np.random.Generator:
        """Get the random number generator used by this object."""

    def max_iter(self, f: np.ndarray, reg_res: rgz.RegularizationResult) -> int:
        """Get the maximum number of iterations for the separation process."""

    def iterative_search(self) -> bool:
        """Determine whether the algorithm should use iterative search."""

Let’s examine each method of the RandomSeparator interface in detail, explaining their roles in the spectral decomposition algorithm and how they interact with each other during execution.

Initial Threshold Calculation#

@abstractmethod
def initial_threshold(self, f: np.ndarray, reg_res: rgz.RegularizationResult) -> float:

Purpose:

This method computes the initial threshold used to identify potential higher-order characters in the first phase of the algorithm. The threshold determines which eigenvalues of the regularized function are considered significant enough to potentially represent higher-order characters.

Parameters:

  • f: The original input function before regularization

  • reg_res: The result of regularizing f, containing eigenvalues and eigenvectors

Returns:

A float value representing the computed threshold

Implementation Notes:

Subclasses should implement this method to return a threshold that balances between:

  • Including all potential higher-order characters (threshold too low).

  • Excluding noise and insignificant components (threshold too high).

Iteration Threshold Calculation#

@abstractmethod
def iteration_threshold(self, f: np.ndarray, reg_res: rgz.RegularizationResult,
                      initial_f: np.ndarray, initial_reg_res: rgz.RegularizationResult) -> float:

Purpose:

This method computes the threshold used during each iteration of the separation process. Unlike the initial threshold, this threshold may adapt to the current state of the algorithm.

Parameters:

  • f: The current function being analyzed.

  • reg_res: The regularization result of f.

  • initial_f: The original input function (for reference).

  • initial_reg_res: The initial regularization result (for reference).

Returns:

A float value representing the computed iteration-specific threshold.

Separation Gap Calculation#

@abstractmethod
def separation_gap(self, f: np.ndarray, reg_res: rgz.RegularizationResult,
                 initial_f: np.ndarray, initial_reg_res: rgz.RegularizationResult) -> float:

Purpose:

This method determines when an eigenvalue is sufficiently separated from the rest of significant eigenvalues. If we denote by \(\Delta\) the output of this method, the function hofa.char.spechoft will consider \(\lambda\) separated if it has multiplicity 1 and there are no other eigenvalues in the interval \((\lambda-\Delta,\lambda+\Delta)\).

Parameters:

Same as iteration_threshold.

Returns:

A float value representing the minimum required gap between eigenvalues.

Implementation Notes:

The separation gap is crucial for determining algorithm termination. Effective implementations consider:

  • The distribution of eigenvalues.

  • The gaps between consecutive eigenvalues.

  • Statistical properties of the eigenvalue spectrum.

Random Number Generator Access#

def get_rng(self) -> np.random.Generator:

Purpose:

This method provides access to the random number generator used internally by the separator. This allows for reproducible results and gives users control over the randomization process.

Returns:

The numpy.random.Generator instance used by this object.

Implementation Notes:

  • The random number generator (RNG) should be properly initialized during object construction.

  • Users can inspect or modify the RNG state if needed (though modifying it may affect reproducibility)

Maximum Iterations#

def max_iter(self, f: np.ndarray, reg_res: rgz.RegularizationResult) -> int:

Purpose:

This method determines the maximum number of iterations allowed for the separation process. It provides a way to control the trade-off between computation time and result quality.

Parameters:

  • f: Current function.

  • reg_res: Current regularization result.

Returns:

The maximum number of iterations to perform.

Implementation Notes:

The default implementation typically returns a fixed value set during initialization. Subclasses can override this to implement dynamic determination based on:

  • The size or complexity of the input data.

  • Computational constraints.

This method ensures the algorithm terminates even if separation is not achieved.

Iterative Search Flag#

def iterative_search(self) -> bool:

Purpose:

This method indicates whether the algorithm should use an iterative approach to find higher-order characters. The iterative approach is generally more efficient than trying to find all components at once.

Returns:

True if iterative search should be used, False otherwise.

Implementation Notes:

Iterative search is typically faster and similarly accurate.

Non-iterative approaches might be preferred when:

  • The number of target components is small.

  • The data has special structure that allows direct decomposition.

This flag is usually set during initialization and remains constant.

Standard Implementation: StandardRandomSeparator#

The StandardRandomSeparator class provides a concrete implementation of the RandomSeparator interface. This implementation is designed to work effectively across a wide range of applications while allowing for customization through various parameters and modes.

Constructor#

def __init__(
    self,
    param_initial: float | int = 1.2,
    mode_initial: str = 'dynamic-relative',
    param_iter: float | int = 1.2,
    mode_iter: str = 'dynamic-relative',
    param_separation: float | int = 0.6,
    mode_separation: str = 'dynamic',
    max_iterations: int = 20,
    iterative_search: bool = True,
    rng: int | np.random.Generator | None = None
):

The constructor initializes the parameters that control the behavior of the separation algorithm across its three main phases: initialization, iteration, and separation. Let us now describe what these options mean.

Initial Threshold Calculation#

def initial_threshold(self, f: np.ndarray, reg_res: rgz.RegularizationResult) -> float:

This method computes the initial threshold for identifying potential higher-order characters based on the input function and its regularization result. The computation uses the information specified during initialization (mode_initial and param_initial).

Mathematical Formulation: The behavior of this method depends on the value of mode_initial, which can take four values. Before presenting them, let us define the following:

  • \(\sigma\) is the standard deviation of the input function f.

  • \(|Z|\) is the size of the base group of f or, equivalently, the product of the function’s dimensions.

The values returned depending on mode_initial are then:

  1. dynamic-relative. This is the default option. First, let \(\rho_{\text{temp}} = \text{param\_initial} \times \sigma \sqrt{\frac{|Z|}{\log|Z|}}\). The value \(\rho_{\text{init}}\) returned is then the middle point of the largest gap of the eigenvalues computed during the initial regularization of \(f\), i.e. in step 1.A of the Minimalistic formulation, which lie in the interval \([\rho_{\text{temp}}/2,\rho_{\text{temp}}]\).

  2. dynamic-strict. The value returned is \(\rho_{\text{init}} = \text{param\_initial} \times \sigma \sqrt{\frac{|Z|}{\log|Z|}}\).

  3. literal-relative. The value \(\rho_{\text{init}}\) returned is the middle point of the largest gap of the eigenvalues computed during the initial regularization of \(f\), i.e. in step 1.A of the Minimalistic formulation, which lie in the interval \([\text{param\_initial}/2,\text{param\_initial}]\).

  4. literal-strict. The value returned is \(\text{param\_initial}\).

Iteration Threshold Calculation#

def iteration_threshold(self, f: np.ndarray, reg_res: rgz.RegularizationResult,
                       initial_f: np.ndarray, initial_reg_res: rgz.RegularizationResult) -> float:

This method computes the threshold used during each iteration of the separation process, following the same mathematical formulation as the initial_threshold but using the parameters specified for the iteration phase (param_iter and mode_iter).

Separation Gap Calculation#

def separation_gap(self, f: np.ndarray, reg_res: rgz.RegularizationResult,
                  initial_f: np.ndarray, initial_reg_res: rgz.RegularizationResult) -> float:

This method calculates the separation gap that determines when eigenvalues are sufficiently separated depending on the value of mode_separation and param_separation.

Mathematical Formulation: Let us use the same conventions as before for \(|Z|\) and \(\sigma\). Then the possible values of mode_separation and the corresponding computation are:

  1. dynamic. This is the default value. In this case the returned value is \( \text{param\_separation} \times \sigma \sqrt{\frac{|Z|}{\log|Z|}}\).

  2. literal. In this case the returned value is \( \text{param\_separation} \).

Random Number Generator Access#

def get_rng(self) -> np.random.Generator:

This method provides access to the internal random number generator, allowing for reproducible results and inspection of the RNG state.

Implementation Notes:

  • The RNG is initialized during construction.

  • Users can provide either a seed or a Generator instance.

  • The same RNG is used for all randomized operations within the class.

Maximum Iterations#

def max_iter(self, f: np.ndarray, reg_res: rgz.RegularizationResult) -> int:

This method returns the fixed maximum number of iterations set during initialization. The input parameters are accepted for consistency with the abstract base class but are not used in this implementation. The returned default value is 20 iterations.

Note: Subclasses can override this method to implement dynamic determination of maximum iterations based on the input data and regularization results.

Iterative Search Flag#

def iterative_search(self) -> bool:

This method returns the value of the do_iterative_search flag set during initialization, which determines whether the algorithm should use an iterative approach for separation. The default value is True.

Performance Considerations:

  • Iterative search is generally faster and similarly accurate.

  • Non-iterative approaches might be preferred for small target component counts.

Result Container: RandomizedSeparationResult#

The RandomizedSeparationResult class provides a structured container for the outputs of hofa.char.spechoft. This NamedTuple encapsulates all the essential information returned by functions that perform spectral decomposition and separation of higher-order characters.

class RandomizedSeparationResult(NamedTuple):
    higher_order_char: np.ndarray,
    separated_eigenvalues: np.ndarray,
    are_separated: bool,
    separation_threshold: float,
    total_iterations: int

Fields#

higher_order_char:

  • Type: numpy.ndarray.

  • Description: The computed higher-order characters after separation.

  • Structure: A multi-dimensional array where the last dimension corresponds to the separated characters.

  • Shape: (f.shape, k) where f.shape is the shape of the original input function and k is the number of separated characters.

separated_eigenvalues:

  • Type: numpy.ndarray.

  • Description: The eigenvalues corresponding to the separated higher-order characters.

  • Structure: 1D array containing the eigenvalues, where the ith eigenvalue corresponds to the eigenvector higher_order_char[...,i]. If there was no iterative search, then this array is in descending order with all eigenvalues separated. Otherwise, they may not be in descending order or separated.

are_separated:

  • Type: bool.

  • Description: Flag indicating whether the algorithm succeeded.

  • Usage: Should be checked before using the results.

separation_threshold:

  • Type: float.

  • Description: The threshold value used to determine successful separation if there was no iterative search or the last threshold used in case of iterative search.

total_iterations:

  • Type: int.

  • Description: Total number of iterations performed by the algorithm.

To conclude, let us demonstrate the usefulness of hofa.char.spechoft with some more complicated examples.

Example: 3 quadratic waves and cubic noise#

Problem statement#

Let \(n=1701\) and consider the following function \(g:\{0,\ldots,n-1\}\to \mathbb{C}\) given by

\[ g(x) = (\exp(2\pi i x^2/n)+\exp(10\pi i t_j x^2/n)+\exp(14\pi i t_j x^2/n))/3. \]

We would like to get the higher-order components (characters) of \(g\). But as usual, imagine that, for \(t=0.1\), instead of \(g\) we are presented with a noisy version \(f(x)=(1-t)g(x)+tf_r(x)\) where \(f_r(x)= \exp(2\pi i x^3/n)\). Note that, even though \(f_r\) is also a polynomial, it has a higher degree, and thus, it should be regarded as noise for the purposes of quadratic regularization.

n = 1701
t = 0.1

x = np.arange(n)
g = np.exp(2 * np.pi * 1j * x**2 / n) + np.exp(2 * np.pi * 5 * 1j * x**2 / n) + np.exp(2 * np.pi * 7 * 1j * x**2 / n) 

f = np.exp(2 * np.pi * 1j * x**3 / n)*t+(1-t)*g/3

Objective#

Try to recover the individual components that make up \(g\).

Strategy#

Simply apply hofa.char.spechoft.

rnd_sep_result = char.spechoft(f, rng=28) 

eichars = rnd_sep_result.higher_order_char

import itertools
import numpy as np
import matplotlib.pyplot as plt

number_of_first_elems_to_plot = 500

v_1 = eichars[:, -1]
coeff_proj_1 = np.mean(f * v_1.conjugate())

v_2 = eichars[:, -2]
coeff_proj_2 = np.mean(f * v_2.conjugate())

v_3 = eichars[:, -3]
coeff_proj_3 = np.mean(f * v_3.conjugate())

# Restrict to plotting window
x_plot = x[:number_of_first_elems_to_plot]

# Target ("green") complex curves
target_1_complex = (1 - t) * np.exp(2 * np.pi * 1j * x_plot**2 / n) / 3
target_2_complex = (1 - t) * np.exp(2 * np.pi * 5 * 1j * x_plot**2 / n) / 3
target_3_complex = (1 - t) * np.exp(2 * np.pi * 7 * 1j * x_plot**2 / n) / 3

targets_complex = [target_1_complex, target_2_complex, target_3_complex]

# Recovered ("blue") complex curves
blue_1_complex = coeff_proj_1 * v_1[:number_of_first_elems_to_plot]
blue_2_complex = coeff_proj_2 * v_2[:number_of_first_elems_to_plot]
blue_3_complex = coeff_proj_3 * v_3[:number_of_first_elems_to_plot]

blues_complex = [blue_1_complex, blue_2_complex, blue_3_complex]

# --- L^2 (squared) error comparison over all 3! possible assignments ---
# Match using the full complex-valued curves

best_perm = None
best_err = np.inf

for perm in itertools.permutations(range(3)):
    total_err = sum(
        np.sum(np.abs(blues_complex[perm[i]] - targets_complex[i])**2)
        for i in range(3)
    )
    if total_err < best_err:
        best_err = total_err
        best_perm = perm

# Reorder blues according to best assignment
matched_blues_complex = [blues_complex[best_perm[i]] for i in range(3)]

# Labels
target_labels_real = [
    r'$\mathrm{Re}\!\left((1-t)\exp(2\pi i x^2/n)/3\right)$',
    r'$\mathrm{Re}\!\left((1-t)\exp(2\pi i\,5x^2/n)/3\right)$',
    r'$\mathrm{Re}\!\left((1-t)\exp(2\pi i\,7x^2/n)/3\right)$',
]

target_labels_imag = [
    r'$\mathrm{Im}\!\left((1-t)\exp(2\pi i x^2/n)/3\right)$',
    r'$\mathrm{Im}\!\left((1-t)\exp(2\pi i\,5x^2/n)/3\right)$',
    r'$\mathrm{Im}\!\left((1-t)\exp(2\pi i\,7x^2/n)/3\right)$',
]

blue_labels_real = [
    r"real part of matched quadratic Fourier character 1",
    r"real part of matched quadratic Fourier character 2",
    r"real part of matched quadratic Fourier character 3",
]

blue_labels_imag = [
    r"imaginary part of matched quadratic Fourier character 1",
    r"imaginary part of matched quadratic Fourier character 2",
    r"imaginary part of matched quadratic Fourier character 3",
]

# --- Plotting ---

# 4 rows x 2 columns:
# top row spans both columns,
# then 3 rows of 2 plots each
fig = plt.figure(figsize=(15, 16))
gs = fig.add_gridspec(4, 2, height_ratios=[1, 1, 1, 1], width_ratios=[1, 1])

# Top plot
ax_top = fig.add_subplot(gs[0, :])
ax_top.plot(
    x_plot,
    np.real(f[:number_of_first_elems_to_plot]),
    linestyle="-",
    color="red",
    marker="o",
    markersize=5,
    linewidth=1,
    alpha=0.6,
    label=r"Real part of original $f$",
)
ax_top.plot(
    x_plot,
    np.real((1 - t) * g[:number_of_first_elems_to_plot] / 3),
    linestyle="-",
    color="green",
    marker="o",
    markersize=5,
    linewidth=1,
    alpha=0.6,
    label=r"Real part of quadratic truth $g$",
)
ax_top.set_xlabel("x")
ax_top.set_title(f"t (Amount of noise): {t:.2f}, n (total length): {n}")
ax_top.set_ylim(-1, 1)
ax_top.legend(loc="upper center", bbox_to_anchor=(0.5, -0.15), ncol=2, frameon=False)

# Bottom 3x2 comparison plots
for i in range(3):
    target = targets_complex[i]
    blue = matched_blues_complex[i]

    # Left: real parts
    ax_real = fig.add_subplot(gs[i + 1, 0])
    ax_real.plot(
        x_plot,
        np.real(target),
        linestyle="-",
        color="lightgreen",
        marker="o",
        markersize=5,
        linewidth=1,
        alpha=0.6,
        label=target_labels_real[i],
    )
    ax_real.plot(
        x_plot,
        np.real(blue),
        linestyle="-",
        color="blue",
        marker="o",
        markersize=3,
        linewidth=1,
        alpha=0.6,
        label=blue_labels_real[i],
    )
    ax_real.set_xlabel("x")
    ax_real.set_ylim(-1, 1)
    ax_real.legend(loc="upper center", bbox_to_anchor=(0.5, -0.15), ncol=1, frameon=False)

    # Right: imaginary parts
    ax_imag = fig.add_subplot(gs[i + 1, 1])
    ax_imag.plot(
        x_plot,
        np.imag(target),
        linestyle="-",
        color="lightgreen",
        marker="o",
        markersize=5,
        linewidth=1,
        alpha=0.6,
        label=target_labels_imag[i],
    )
    ax_imag.plot(
        x_plot,
        np.imag(blue),
        linestyle="-",
        color="blue",
        marker="o",
        markersize=3,
        linewidth=1,
        alpha=0.6,
        label=blue_labels_imag[i],
    )
    ax_imag.set_xlabel("x")
    ax_imag.set_ylim(-1, 1)
    ax_imag.legend(loc="upper center", bbox_to_anchor=(0.5, -0.15), ncol=1, frameon=False)

fig.tight_layout()
plt.show()
../_images/20d9ecadee23794ab114b3beff8806d12c64b39fa23bd1d4403224193116feb3.png

Example: Locating quadratic behaviour in data analysis#

The next example shows how hofa.char.spechoft can locate quadratic behaviour spatially, not only if it is given as a superposition of quadratic (or higher-order) waves.

Problem statement#

Let \(n=1701\), \(m=500\), and let us define \(f:\{0,\ldots,n-1\}\to \mathbb{C}\) as follows:

\[\begin{split} f(x) = \begin{cases} \exp(2\pi i x^2/n) & \text{if }x\le m\\ \exp(10\pi i x^2/n) & \text{if }x> m. \end{cases} \end{split}\]
n = 1701
m = 500

x = np.arange(n)
f = np.exp(2 * np.pi * 1j * x**2 / n)
f[m:] = np.exp(2 * np.pi * 5 * 1j * x**2 / n)[m:]

Strategy#

We would like to isolate as much as we can the two components of this function.

To do so, we apply hofa.char.spechoft. This time, for illustration purposes, we add a callback function to print the advances in the iterations.

rnd_sep_result = char.spechoft(f , rng=28, callback = char.PrintRandomizedSearchCallback())

eichars = rnd_sep_result.higher_order_char

Total number of iterations 1, iterative search = True, eig vals found = 2, target num eig vals = 2

Results#

We can already see from the message above that hofa.char.spechoft was able to correctly discover that there are two markedly different components hidding in \(f\). First, let us print the absolute value of the higher-order characters found. Note that it clearly separates the two quadratic components.

fig, ax = plt.subplots(figsize=(12, 5))

colors = ['blue', 'green', 'yellow', 'red']

for i in range(2):
    # Plot recovered function in blue
    v = rnd_sep_result.higher_order_char[...,i]
    coeff_proj = np.mean(f*v.conjugate())
    ax.plot(x, np.abs(coeff_proj*v), linestyle='-', color=colors[i], marker='o', markersize=3, linewidth=1, alpha=0.6, label='abs value of eigencharacter')

# Customize the plot
ax.set_xlabel("x")
#ax.set_ylim(-1,1)

# Legend below the plot
fig.subplots_adjust(bottom=0.25)  # Adjust layout to fit legend below
ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.15), ncol=2, frameon=False)

# Show the plot
plt.show()
../_images/825ba56ccfe84316d75141552e9c04ec9803ee63d94198ae26ab1806b6778bd8.png

Let us plot below the different components together with the ones found by the spectral higher-order Fourier transform.

import itertools
import numpy as np
import matplotlib.pyplot as plt

number_of_first_elems_to_plot = 1000

v_1 = eichars[:, -1]
coeff_proj_1 = np.mean(f * v_1.conjugate())

v_2 = eichars[:, -2]
coeff_proj_2 = np.mean(f * v_2.conjugate())

# Restrict to plotting window
x_plot = x[:number_of_first_elems_to_plot]

# Target ("green") complex curves
target_1_complex = np.zeros_like(x_plot, dtype=complex)
target_1_complex[:m] = np.exp(2 * np.pi * 1j * x_plot**2 / n)[:m]

target_2_complex = np.zeros_like(x_plot, dtype=complex)
target_2_complex[m:] = np.exp(10 * np.pi * 1j * x_plot**2 / n)[m:]

targets_complex = [target_1_complex, target_2_complex]

# Recovered ("blue") complex curves
blue_1_complex = coeff_proj_1 * v_1[:number_of_first_elems_to_plot]
blue_2_complex = coeff_proj_2 * v_2[:number_of_first_elems_to_plot]

blues_complex = [blue_1_complex, blue_2_complex]

# --- L^2 (squared) error comparison over all possible assignments ---

best_perm = None
best_err = np.inf

for perm in itertools.permutations(range(2)):
    total_err = sum(
        np.sum(np.abs(blues_complex[perm[i]] - targets_complex[i])**2)
        for i in range(2)
    )
    if total_err < best_err:
        best_err = total_err
        best_perm = perm

# Reorder blues according to best assignment
matched_blues_complex = [blues_complex[best_perm[i]] for i in range(2)]

# Labels
target_labels_real = [
    r'$\mathrm{Re}\!\left(\exp(2\pi i x^2/n)/2i\right)$ for $x< m$',
    r'$\mathrm{Re}\!\left(\exp(10\pi ix^2/n)/2i\right)$ for $x> m$',
]

target_labels_imag = [
    r'$\mathrm{Im}\!\left(\exp(2\pi i x^2/n)/2i\right)$ for $x< m$',
    r'$\mathrm{Im}\!\left(\exp(10\pi ix^2/n)/2i\right)$ for $x> m$',
]

blue_labels_real = [
    r"real part of matched quadratic Fourier character 1",
    r"real part of matched quadratic Fourier character 2",
]

blue_labels_imag = [
    r"imaginary part of matched quadratic Fourier character 1",
    r"imaginary part of matched quadratic Fourier character 2",
]

# --- Plotting ---

fig = plt.figure(figsize=(15, 16))
gs = fig.add_gridspec(4, 2, height_ratios=[1, 1, 1, 1], width_ratios=[1, 1])

# Top plot
ax_top = fig.add_subplot(gs[0, :])
ax_top.plot(
    x_plot,
    np.real(f[:number_of_first_elems_to_plot]),
    linestyle="-",
    color="red",
    marker="o",
    markersize=5,
    linewidth=1,
    alpha=0.6,
    label=r"Real part of input $f$",
)
ax_top.set_xlabel("x")
ax_top.legend(loc="upper center", bbox_to_anchor=(0.5, -0.15), ncol=2, frameon=False)

# Bottom 3x2 comparison plots
for i in range(2):
    target = targets_complex[i]
    blue = matched_blues_complex[i]

    # Left: real parts
    ax_real = fig.add_subplot(gs[i + 1, 0])
    ax_real.plot(
        x_plot,
        np.real(target),
        linestyle="-",
        color="lightgreen",
        marker="o",
        markersize=5,
        linewidth=1,
        alpha=0.6,
        label=target_labels_real[i],
    )
    ax_real.plot(
        x_plot,
        np.real(blue),
        linestyle="-",
        color="blue",
        marker="o",
        markersize=3,
        linewidth=1,
        alpha=0.6,
        label=blue_labels_real[i],
    )
    ax_real.set_xlabel("x")
    ax_real.legend(loc="upper center", bbox_to_anchor=(0.5, -0.15), ncol=1, frameon=False)

    # Right: imaginary parts
    ax_imag = fig.add_subplot(gs[i + 1, 1])
    ax_imag.plot(
        x_plot,
        np.imag(target),
        linestyle="-",
        color="lightgreen",
        marker="o",
        markersize=5,
        linewidth=1,
        alpha=0.6,
        label=target_labels_imag[i],
    )
    ax_imag.plot(
        x_plot,
        np.imag(blue),
        linestyle="-",
        color="blue",
        marker="o",
        markersize=3,
        linewidth=1,
        alpha=0.6,
        label=blue_labels_imag[i],
    )
    ax_imag.set_xlabel("x")
    ax_imag.legend(loc="upper center", bbox_to_anchor=(0.5, -0.15), ncol=1, frameon=False)

fig.tight_layout()
plt.show()
../_images/d6d8c4a900f145c0e4134b3d4cb448b265e47c77169b9afa1245a41d7f85a0a9.png
On this page

Table of Contents

Show Source