As a symbolic paradigm, many HDC algorithms operate on a set of dissimilar atomic hypervectors that are randomly generated and near-orthogonal, assuming that symbols are not related at all. However, the assumption will not always be appropriate in practical tasks. Therefore, we have seen HDC algorithms, when handling bio-signals and images, explicitly manipulate the similarity among atomic hypervectors such as maintaining a discrete set of similarity levels. However, the manual assignment of these hypervectors can be problematic, which inevitably causes information loss during quantization, not to mention that such an arbitrarily assumed similarity relationship may not be helpful for learning.

To explain how the VFA encoding captures the relation between data, we note that this encoding coincides with the Random Fourier Features (RFF) encoding, an efficient approximation of kernel methods. The following theorem by Salomon Bochner (1899–1982) serves as the foundation for this well-defined similarity relationship (Rudin, 2017).

Theorem 1 (Bochner). For any continuous shift-invariant and positive definite kernel K(x1−x2):ℝM→ℝ, there must exist a non-negative measure p(ω) such that K is the Fourier transform of a non-negative measure p(ω). Additionally, if K is properly scaled, p(ω) is a proper probability measure.

The proof of this theorem is provided in Rudin (2017). If we assume ζ(x) = e^jωTx, then Theorem 1 leads to the following equation: K(x1−x2)=∫ℝMp(ω)ejωT(x1−x2)=Eω[ζ(x1)ζ(x2)¯]. This means that, with the correspondence between kernel K and measure p(x), we can transform original inputs to a space where dot products are unbiased estimates of kernel similarities. In other words, there exist sequences of transformations ϕD:ℝM→CD such that  ϕD(x1)TϕD(x2) converges uniformly to the given kernel K(x₁−x₂):

Rahimi and Recht (2007) proposed an alternative set of Random Fourier Features (RFF) such that the components of the encoded vectors are real and the kernel approximation converges equally fast. To construct a real-valued RFF, we can leverage this high-dimensional mapping for HDC encoder as the following:

where cos. represents the element-wise cosine function, Ω∈ℝ^D×M is a randomly generated encoding matrix, and b∈ℝ^D is the offset hypervector. Row vectors in Ω are generated by drawing D i.i.d. samples ω₁, …, ω_D from p(ω) and elements in b are sampled from U[0,2π]. The random distribution p is selected through Theorem 1 given a preferred kernel K(Δ). In other words, the probability density is calculated with a Fourier transform: p(ω)=12π∫ℝMexp(-iωTΔ)K(Δ)dΔ. For example, previous HDC work (Hernández-Cano et al., 2021) uses Normal distribution N(0,1) for p(ω) since it wants to approximate the Gaussian RBF kernel.

There are key lessons from the theory of VFA encoding discussed above:

While the current VFA method brings many benefits, the biggest drawback of this method is that ϕ_D(x) is essentially a static mapping, which makes the encoding less adaptive. Our work aims to leverage the insight from Bochner's theorem to learn the kernel adaptively through its random Fourier features.

Ideally, we expect the HDC encoder to provide a useful high-dimensional representation that helps separate the data points for classification or linearize the inherent non-linear regression tasks. Particularly in hyperdimensional regression, we are interested in finding the best hypervector w∈ℝ^D such that the linear regression after encoding y(x) = ϕ(x)^Tw are, in average across the training set, as close as possible to the true labels in terms of ℓ₂ norm. Additionally, we introduce an ℓ₂ regularization coefficient λ to get more stable estimators. Thus the loss function is the dampened least squares, which can be expressed as

where y=(y1,y2,…,yN)N are the known response variables, and Z = Φ(X; Ω, b)∈ℝ^N×D are the encoded hypervectors of the input data x₁, …, x_N:

In previous approaches for HDC regression (Hernández-Cano et al., 2021; Kleyko et al., 2021), the model hypervector w is learned in an iterative fashion, where hypervectors are bundled together guided by regression errors. However, they face issues with proper hyperparameter selection to achieve the highest prediction quality. On the other hand, the loss function in Equation (3) has a known minimizer w^ which is known as the ridge estimator:

In this paper, we leverage this statistical approach to obtain better stability during learning and more direct parameter tuning.

As we mentioned in the previous section, the encoding in VFA, the closed solution in Equation (5) has an assumption that the regression problem on Z is linearly solvable, as the result of mapping to hyperspace. However, for an arbitrary regression task, it is very likely that the static VFA encoding (due to the fixed K(Δ) and p(ω)) becomes sub-optimal. This work looks to address this problem by presenting an adaptive HDC encoder design.

Care must be taken when designing the HDC encoder, as it needs to ensure that the appealing properties of the hypervectors are sustained. To ensure the holographic representation, we require the randomized instantiation of the encoder, and thus we cannot directly perform gradient descent upon an instantiated encoding matrix, as it may destroy both: information about the data may be distributed locally and partially in the output vector, and the trained encoder have minimum randomization.

To circumvent this, we take inspiration from the VFA work that highlights the encoder-kernel correspondence and the importance of selecting a proper distribution p(ω). Recall from Section 3.1 that there is a correspondence between kernel K and measure p(x). In addition, it induces families of encoders {_ϕD}D∈ℕ parameterized by the samples from the distribution such that the inner product between the encoded vectors approximates the kernel. This approximation improves increasingly well with the dimension of the encoder (codomain). By learning the distribution from which we sample the encoding matrix, we will be able to construct an adaptive HDC encoder that provides a more suitable hyperdimensional representation, adding to existing appealing HDC properties.

In FLASH, we define a parametric family of functions F={fθ:fθ:ℝM→ℝM,θ∈Θ}, which will be referred to as generators. Upon receiving random inputs, it can be used to sample random vectors ω₁, …, ω_D required in the encoding function (Equation 4). Compared to parameterized distributions such as the Gaussian reparameterization approach, this approach encapsulates a more expressive family of distributions. Because of Bochner's theorem, exploring F is equivalent to exploring a corresponding set of continuous shift-invariant kernels, and thus the encoding family we consider is expected to be rich as long as F is.

Inherited from the RFF methods, a key benefit of this arrangement is that FLASH encoding does not require us to use explicitly the kernel K(Δ) function, nor the probability density p(ω), nor to perform expensive Fourier transforms. Our goal is to find f∈F that, with a high probability, gives the optimal or near-optimal encoding matrix of solving the regression problem at hand; this ensures the quality and robustness of the encoder that it generates.

To learn the distribution efficiently, we restrict F to be a family of fθ:ℝM→ℝM differentiable neural networks, i.e., the network input size equals the output size. To sample ωs using f_θ, we first draw a random vector ϵ~N(0,I) as the input, and then obtain a transformed random vector using ω = f_θ(ϵ). Sampling D i.i.d. random noises and passing them through the generator can give us a matrix of base hypervectors (i.e., Ω), which can be used to perform the encoding. This can be understood as a generalized reparameterization, where we learn a surrogate function f_θ(ϵ) instead. Note that we chose to sample noise vectors from the standard normal distribution for convenience, but different choices can be made as well. This architecture gives us a very rich family of generators F, which are cheap to evaluate and cheap to train, as we will see in the next few sections. In addition, because ϵ is a random vector, ω = f_θ(ϵ) is one too, which means that there exists a probability density (or mass) function p_θ(ω) for each generator f_θ.

In FLASH, we aim to maximize encoder performance in generating a good HDC representation of our data for the regression task. In particular, we opt to evaluate the learned encoder using the loss function proposed in Equation (3). However, because random sampling is involved at the moment of generating the encoding, we chose to minimize the expected value instead. Note that this expected value is taken with respect to all the possible encoding Φ, whose encoding matrix Ω and offset hypervectors b are randomly sampled. Thus, in Equation 6, we seek to find the parameters in f_θ that minimize the following loss term for adapting the HDC encoder:

where LR(w) is the regression loss defined in Equation (3). In Equation 7, using the ridge estimator w^ and the law of the unconscious statistician, we can expand the previous expectation term to:

where E is the matrix containing the D random vectors ϵ₁, …, ϵ_D. Thus, we can obtain an unbiased estimator of the encoding loss using a simple Monte Carlo estimator with a single sample. That is, if we sample E and b from their respective distributions, we obtain an unbiased estimator of LE(θ):

Note that it is possible to sample multiple noise matrices E₁, E₂, … to lower the variance of the estimator, but in order to accelerate the computation we don't explore this alternative.

We explore the parameter space of the generator f_θ using stochastic gradient descent, where θ←θ-η∇θL^E(θ). Its easy to show that ∇θL^(θ) is an unbiased estimator of the true gradient ∇θL(θ) because the expectation in LE does not depend on the parameters θ. It is important to note that L^E is differentiable with respect to its parameters θ. Indeed, every operation shown in Equation (8) is well behaved: Ω = f_θ(E) is clearly differentiable, and so is Z = Φ(X; Ω, b) because φ(x;Ω,b)=2Dcos.(Ωx+b).

Until here, we have covered how to parameterize and learn the sampling distribution. This equips FLASH with an encoding module well-adapted. Still, people may wonder if learning the encoding matrix Ω and offset b through gradient descent is a good alternative. We acknowledge that this might be a more direct measure, however, it will jeopardize the holographic property of HDC since it cannot guarantee that the encoding matrix is i.i.d. This means that the information in the encoded hypervector is no longer evenly distributed, and errors or noise in the encoding process will lead to higher performance loss due to the lack of hyperdimensional redundancy. Our proposed measure will ensure that FLASH will maintain the holographic HDC representation after tuning the encoder.

The training in FLASH is a two-stage process where we first learn the generator f_θ(ϵ), i.e., in place of sampling from p_θ(ω)). Based on the first training stage, we then perform the model training that gives the regression hypervector w. In the second stage, the encoder will be generated using f_θ and remain static as it has been optimized. Notice that the dimensionality D can be the same or different in these two training stages. As mentioned in the previous section, our optimization method for generator training is well-defined; but in practice, it can be further approximated to obtain faster convergence. Several algorithms have been widely used to accelerate the computation of ridge estimator (Paige and Saunders, 1982; Defazio et al., 2014). However, computing LR(θ) at every iteration for encoder learning ends up adding up the overhead. Recall that in prior HDC algorithms (Hernandez-Cane et al., 2021; Hernández-Cano et al., 2021; Ni et al., 2022a), D is supposed to be a high dimensionality such that model hypervectors have a larger capacity. In FLASH, we instead propose to decouple the high dimensionality requirement from the encoder/generator training since the generator f_θ itself operates in ℝ^M. When training f_θ, we encode data to ℝ^D′ instead, where D′ < D in order to accelerate the process. As for the regression process, we use a slightly larger dimensionality D for better regression accuracy. In fact, adapting the HDC encoder at first will also lower the requirement for model hypervector dimensionality D and thus reduce the training cost. Our results in the experiment show that FLASH, with a lower dimensionality, has comparable regression quality to the prior HDC method.

In this section, we discuss the time complexity of training our proposed method. Below we describe at a high level the steps required in our approach.

In the first step, the overhead of computing LE is considered minor as we encode data to D′-dimensional space, limiting the cost of computing the estimator w^. Generating random bases ω_i = f_θ(ε_i) requires D forward passes of f_θ, which will give a hyperdimensional mapping.

Encoding the data Z=Φ(x;Ω,b)=2Dcos.(xΩT+b), requires O(NMD) operations, corresponding to the asymptotic of the most taxing operation - matrix multiplication of the N×M matrix X with a M×D matrix Ω^T, where N is the number of samples and M the number of features in original space.

The last step, computing w = (Z^TZ+λI)−1Z^Ty is, according to experimental evaluation, the most time-consuming step in our design. The theoretical time complexity is dominated by the Z^TZ multiplication and the inverse operation, requiring O(ND2) and O(D3) operations, respectively. Thankfully, the dampened linear least squares loss, Equation 3), has been heavily studied and several algorithms that approximate its solution exist such as LSQR (Paige and Saunders, 1982). Moreover, the experimental evaluation suggests that with relatively small values of D (500 ≤ D ≤ 2000) we can obtain very accurate predictions (as shown in Figure 5); with this configuration, we observe linear training time in the number of samples.

In this section, we show that minimizing LE with respect to θ is equivalent to maximizing a lower bound of the log-likelihood of the posterior distribution with joint parameters (θ, w).

Recall that given a dataset of independent observations D, random parameters θ∈Θ with prior distribution θ~p(θ) the posterior distribution is p(θ∣D), which by Bayes' Theorem can be expressed as p(θ∣D)∝p(D∣θ)p(θ), where p(D∣θ) is called the likelihood.

In regression analysis, we assume the relation y = Xw+ε, where ε is a random unobserved noise. Ordinary least squares regression sets the likelihood to be normally distributed y∣w~N(xw,σ2I). As shown in Equation 9, maximizing the log-likelihood is equivalent to minimizing the squared error loss:

In order to work in the high-dimensional space, we must add the encoding to the equation, or in this case the generator parameters θ. We instead assume the relation y = Zw+ε, where Z = Φ(X; Ω, b). Because the regression coefficients depend on Ω and b, we use the conditional likelihood y∣Ω,b,w~N(Zw,σ2I). We add the independence assumption between θ and b. In Equation 10, the maximum log posterior is shown to be bounded using Jensen's inequality:

We implement the proposed design using Python on the Intel Core i7-12700K CPU platform. The core process of adapting the encoder is implemented using Pytorch, and the regression process is based on the implementation provided by Scikit-Learn. We evaluate the accuracy of our design on several practical regression datasets listed in Table 1, with up to 20,000 samples and 80 features. Table 2 describes the baseline models used to compare with our design, including ridge regression that also leverages RFF approximation and the previous state-of-the-art HDC-based regression algorithm RegHD. During the experiments, we test two settings for our design (FLASH and A-FLASH), slightly different in model size and dimensionality. The name of the second setting stands for “Accurate FLASH,” which has larger dimensionality and model size. In Table 3, we provide the hyperparameters used for the different regression models in our experiments.

In this section, we analyze the performance of the proposed design in custom 1D regression problems of the form y = f(x)+ε with ε~N(0,1) and different choices of target function f. In Figure 2, we also show the encoder's probability distribution p_θ(ω) learnt in the process and the actual kernel function K(Δ) = E[ϕ(x)^Tϕ(0)]. For comparison, RBF kernel has associated a Gaussian N(0,1γI) distribution.

From this experiment, we conclude that our optimization proposal for the encoder loss LE works well in practice and the shapes of learned distributions are varied for each dataset. Moreover, we observe that our proposal adapts to different scales in the data making a clear distinction with SVR. For instance, the first predicted function f(x) = 10x², where SVR clearly underperforms where |x|>1.1.

In this section, we compare the performance of FLASH (as well as A-FLASH) against several baseline regression algorithms using multiple regression datasets. We perform 5 times repeated 5-fold cross-validation in each dataset and report the average prediction quality, confidence intervals, statistical tests for significance, and runtime taken for each fold. We select the most important hyperparameters in SVR and our design using grid search. Our results are summarized in Figure 3.

We observe that our approach is always comparable in accuracy with other state-of-the-art approaches. The accurate version of our approach (A-FLASH) is consistently ranked at the top. Particularly, because our encoder is learnable and well-adapted, we are generally more accurate than other algorithms leveraging static encoder or fixed kernel. In comparison with the prior HDC-based method, A-FLASH achieves significantly better quality without adding notable overhead for inference. In addition, the fast version of our approach (FLASH) is generally among the fastest models. During inference, it is faster than other baselines, including classical kernel-based approaches such as SVR and Kernel Ridge. This is because our prediction complexity is constant with respect to the number of samples. On average, FLASH is about 3.7 × faster inference than the RegHD, 5.5 × faster than kernel ridge/RFF ridge, and 13.75 × faster than SVR.

In this section, we create the Friedman regression datasets (Friedman, 1991) with an increasing number of samples to test the scalability of the proposed algorithm and compare it with other approaches. We observe that our approach is well-suited for large-scale data as we have a linear trend in the time taken to train and also inference time. Meanwhile, the time taken to train classical kernel approaches such as SVR and kernel ridge grows noticeably faster due to their higher computational complexity. Our results are summarized in Figure 4. The leftmost plot shows that our approach is the fastest to achieve high prediction quality even with a small number of samples; in fact, FLASH constantly achieves better accuracy when the training set grows. In terms of inference speed, FLASH is about twice as fast as the other approaches with 5000 training samples, and the gap in between continues to expand.

In this section, we explore the impact of dimensionality (D) in our design for various datasets. Figure 5 displays our results in terms of prediction quality (MSE) and time taken to train the model for different values of D. In the section on “Time Complexity,” we derived the time complexity of our approach to be O(ND2), which is consistent with our experimental results. However, it is worth mentioning that even for relatively small dimensionality (e.g., D = 500) the gain of accuracy for further increasing dimensionality is not significant. Thus, even if the theoretical complexity of the approach is large, in practice, we can obtain acceptable results rapidly.