hyperdimensional-computing / torchhd Goto Github PK

We can use the same design PyTorch uses, we can extend the dtypes, and extend the different tensor types i.e. FloatTensor. Then in the HDC operations we can check the instance type of the hypervector to change the behavior.

import torchhd

hv = torchhd.functional.random_hv(10, 1000)  # bipolar by default (torch.float)
hv = torchhd.functional.random_hv(10, 1000, dtype=torch.bool)
hv = torchhd.functional.random_hv(10, 1000, dtype=torch.complex64)

torchhd.functional.bind(hv[0], hv[1])  # works for any datatype

Given the existence of alternate versions of the primitive functions, e.g. binding by circular convolution or elementwise product, (see #73), people will want to compare them.

Ideally, a user would write something like: similarity(bundle(bind(A,B),C), X)
and be able to execute it for different definitions of the functions without having to rewrite that code snippet.

The implementation might be as simple as something like
similarity = similarity.cosine; bind = bind.convolution; bundle = bundle.add
before executing the code snippet. (Sorry if that's not pythonic.)

At the moment {Torchhd} implements bind(), bundle(), etc. which sort of implies that the current definitions are the only definitions. I think the names bind, bundle, etc. should be reserved for the generic operators and the actual implementations should have more specific names like bind.convolution and bundle.add. Then the specific operators can be assigned to the generic names as needed (in a pythonic way).

I think this approach should be taken to all the functions that might be viewed generically: definitely VSA primitive operators and similarity; very probably the encodings and atomic vector generators; probably cleanup.

Maybe this is absolutely trivial in python and doesn't require any design as such. If that's the case, then the documentation should explicitly show how to compare alternate implementations in a pythonic way.

Explicitly state that we are happy to have researchers add their new algorithms, and encoding schemes

generator = torch.Generator()
generator.manual_seed(seed)
hv = functional.random_hv(10, 10000, generator=generator)
G = structures.Graph(dim_or_input=10000, directed=True)

    G.add_edge(hv[0], hv[1])
    G.add_edge(hv[0], hv[2])
    G.add_edge(hv[1], hv[2])

    print(functional.cosine_similarity(G.node_neighbors(hv[1]), hv))

tensor([-0.0142, 0.0021, -0.0057, -0.0036, -0.0141, -0.0020, 0.0037, -0.0062,
0.0059, -0.0022])

We would like to implement as part of this library the train, test, and predict utility functions like those provided in: https://github.com/rishikanthc/torch-hd/blob/main/torch_hd/utils.py. This issue serves as the discussion thread for how to design theses functions to be simple, and generally applicable for HDC learning applications.

Add initialization methods to create graphs, and trees, etc. from a fully provided definition instead of one-by-one.

• Each class should also be initializable with a single hypervector, i.e., Sequence(Tensor) where the dimensions, dtype, and device are then inferred from the Tensor.
• Sequence.from_tensor(Tensor)
• Graph.from_edges(Tensor) the tensor should follow the shape of the edge hypervectors used in Torch Geometric

Since the embeddings are designed to be compatible with torch.nn modules it makes sense to give them the same Python module name.

The current implementation of bind is self-inverse (because bipolar and boolean hypervectors are equivalent to complex hypervectors with phase quantised to 2 levels). In general bind is not self-inverse, e.g. in HRR, bind is implemented as circular convolution and unbind as circular correlation.

However, rather than implement unbind as the inverse of bind I think it's more helpful to have a unary "multiplicative"-inverse operator that applies to hypervectors. Then unbinding is implemented as binding with the multiplicative inverse of one of the arguments.

bind(inv_mul(A), A) = 1 (where 1 is the binding identity element)

This makes the algebraic nature of VSAs clearer than using an unbind operator.

Likewise, you should also implement a unary "additive"-inverse operator on hypervectors.

bundle(inv_add(A), A) = 0 (where 0 is the bundling identity element, although 0 isn't explicitly represented in every VSA system)

See https://rgayler.github.io/VSA_altitude_hold/vsa_basic_operators.html#3_Negate_vector for an example of the additive inverse.

We will move the discussion of binary hypervector implementation to this issues such that #25 can focus on the complex representation and varying bit-widths.

hv = torchhd.random_hv(10, 1000, dtype=torch.bool)  # creates a boolean tensor
torchhd.bind(hv[0], hv[1])  # performs XOR
torchhd.bundle(hv[0], hv[1], hv[2])  # needs three inputs to calculated the majority function

Add two utility functions that will generate a similarity plot. There first one gets the similarity characteristics of a basis hypervector set. The second shows the similarity of one hypervector with a set.

• Set characteristics
• Hypervector with set similarities

Bind and bundle are usually defined as binary operators with unweighted operands. This reflects the immensely strong mathematical preference for unary and binary operators and the fact that most multi-argument operations can be expressed as a composition of unary and binary operators (the higher-order operators being strictly redundant in this case).

In many cases the generalisation of bundle more than 2 operands is pretty obvious (and the generalisation to weighted operands is possibly less obvious). You can generally achieve the same result by composing multiple binary bundle operators, but this can be less convenient and is possibly more computationally expensive. Bundling operators are generally non-associative. That is, bundle(bundle(A,B),C) != bundle(A,bundle(B,C)). This is easy to see in randsel bundling because of the conversion of weights to sampling probabilities per application of the operator. Thus bundle.randsel(bundle.randsel(A,B),C) ~= bundle.randsel(1*A,1*B,2*C) and bundle.randsel(A,bundle.randsel(B,C)) ~= bundle.randsel(2*A,1*B,1*C), where ~= should be interpreted as meaning preserving equal similarity relations to the operands.

bundle.randsel with k operands only needs to generate one random hypervectors of values from 1 to k to select between the k operands. Composing the equivalent operation from binary operators requires one random hypervector to be generated per binary operator, so is more computationally expensive. Also, if a k-many bundle is natural for the problem you are solving it is mildly tedious to work out the correct weights for the binary operators to get the desired weights for the k-operand operator.

The scalar weights are quite handy for constructing target hypervectors to compar with calculated hypervectors. Also, using weighted bundling is a natural way to create continuous level encodings.

There is probably less call to have a k-operand weighted binding operator because binding operators are generally associative. However, it may be the case in some problems that it is the natural form to take, so why not implement it if it's easy? For example, see https://rgayler.github.io/VSA_altitude_hold/vsa_basic_operators.html#5_Multiply_vectors

We should make a page in the documentation that gives those new to HDC a place to start. For example:

• Unfamiliar with HDC: start by watching Kanerva's lecture
• Familiar with the basics: some simple paper to get them started

Currently we have 3 different similarity functions:

1. hamming_similarity
2. cosine_similarity
3. dot_similarity

And with the future introduction of complex hypervectors we will likely add a forth one if we follow the current design. I think, however, that we should only provide one similarity function that changes it's behavior based on the dtype of the input tensors. It would also be nice if it handles batched operations, i.e., with input shapes (*, d) and (n, d) the output shape should be (*, n) which has the similarity score for each input sample against each other element.

In order to unify the output domain we can stick to the [-1, +1] range that the cosine similarity and the complex variant of cosine similarity produce where 0 means orthogonal, +1 the same, and -1 the exact opposite. We can simply scale the hamming_similarity to fall in this domain.

The dot_similarity will then be removed from the library but is still available as part of PyTorch. And can therefore still be used in specific instances.

API design

x = torchhd.random_hv(10, 10000)
torchhd.functional.similarity(x, x)  # aliased as torchhd.similarity(x, x)

In the replace function, right now passing the old value is needed, we could do the lookup ourselves with the key and avoid having the user giving us the old value

Using torch.prod() we can efficiently implement a binding-based sequence. Moreover, we can perhaps improve the n-grams function to use torch.prod instead of a for loop.

Use cases

Binding-based sequences are used in n-grams, binary-trees, and finite state automata use them, they can also be used for efficient cross-product calculation (see #6), the cookbook uses them to represent tuples but I'm not sure about that use case.

Implementation

• functional.combination(Tensor) binds all the hypervectors on the -2nd dimension
• functional.sequence_permute(Tensor) shifts all hypervectors in tenser by a decreasing amount such that the last one is not permuted, useful reference: https://discuss.pytorch.org/t/how-to-roll-some-rows-of-a-matrix/36722.
• functional.distinct_sequence(Tensor)

It would be good to support TorchScript so HDC models can be exported to be run in a production environment.

Since text data is often used, we can add a wrapper around the Random embeddings to store a mapping of characters to random hypervectors (Text(chars, padding=" ")).

Paper: https://arxiv.org/pdf/2106.05268.pdf
GitHub repo: https://github.com/denkle/HDC-VSA_cookbook_tutorial

For example the functional.ngrams function can be implemented in a number of ways, we should look into which way is the fasted both on CPU and GPU. The same goes for other functions, classes and loading of datasets.

The current implementation keeps track of a Python list of PyTorch tensors representing the hypervectors, we should do some research into whether dynamically allocating a tensor of a variable size such that all hypervectors can live in one PyTorch tensors would improve the performance while not hurting usability.

The current implementation of the permutation operator is shift/rotate. This is only one permutation. There is a shift argument that specifies how far to shift, but all the values of that correspond only to different powers of the one permutation operator.

Some applications require more than one permutation operator. For example, when encoding a ternary tree you might use a different permutation to encode each of the children of each node.

With hypervectors a randomly selected permutation is with exceptionally high probability orthogonal to any other randomly selected permutation. Each permutation is fixed for the life of a VSA circuit, so each random permutation is generated at-initialisation. You can have a generic permutation function that takes a permutation vector (specifying the permutation to be applied) as an argument. For example, see https://rgayler.github.io/VSA_altitude_hold/vsa_basic_operators.html#4_Permute_vector
Any VSA circuit generally uses a small number of fixed distinct permutations.

When you use a permutation you often need its inverse, so a function to create the inverse of a given permutation is generally needed (see https://rgayler.github.io/VSA_altitude_hold/vsa_basic_operators.html#42_Generate_inverse_permutation).

In some applications you need low powers of a given permutation. The permutation vectors for these can be precalculated rather than applying the permutation operator multiple times. Calculating the inverse permutation vector is a special case of calculating an arbitrary power of a permutation - so it might be good to create a generalised function for v\calculating powers of permutation vectors.

This issue is to discuss design improvements around the access to the exact previous hypervector value in the data structures. Right now we require the user to pass the previous version for every mutation method. We can think of designs to provide this behavior. For instance:

hv = torchhd.random_hv(10, 10000)
S = torchhd.structures.Sequence.from_tensor(hv)
S.replace(2, hv[2], hv[5])

Could be:

hv = torchhd.random_hv(10, 10000)
S = torchhd.structures.Sequence.from_tensor(hv)
S.replace(2, hv[5])  # not passing the old value

however this requires the data structure to have access to the exact hypervector. The discussion here is how to implement that in a way that give ample freedom to the user to experiment with various cleanup memories.

This is a cool project and thanks for referencing our paper for random-projection encoding methods! Just FYI - the random projection encoding method you have implemented is a bit different from those discussed in the paper referenced. Let x be an $n$-dimensional input point to encode. I think the particular embedding you are referring to is the following:

z = sign(Mx)

where M is a d x n dimensional matrix whose rows are sampled from the uniform distribution over the unit-sphere. A simple way to generate samples from the uniform distribution over the unit-sphere is to normalize a sample from the $n$-dimensional standard normal distribution (see here for more info). Note that I'm not sure sampling from the uniform distribution over [-1,1]^n and normalizing results in a uniform distribution over the sphere. In particular, I think this approach results in too little mass distributed near the equator and poles. The following code should do the trick:

# Generate the embedding matrix:
import numpy as np
d = 10000; n = 100
M = np.random.normal(size=(d,n))
M /= np.linalg.norm(M, axis=1).reshape(-1,1)

# encode a point:
x = np.random.rand(n)
z = np.sign(M.dot(x))

The sign function is important because the sense in which the encoding preserves distances between points is different without it (and I'm not sure is what one would want). You may not want to use the sign function because it messes with gradients (e.g. its derivative is zero everywhere except at zero, where it does not exist). If you want to omit the sign function and use a linear projection, I would recommend looking into the "Johnson-Lindenstrauss Transform" (see here).

We should try to find a more appropriate term for the binding-based sequence.

Possible alternatives:

• opaque_sequence
• unique_sequence
• diversity_sequence
• contrasting_sequence
• bindlist with sequence as bundlelist or multisetlist

Allow for multiple, different versions of the primitive VSA operators (bind, bundle, permute, and maybe others).

This is needed, for example to introduce HRR binding, and randsel bundling (#75).

This is related to #72 and #25, which appear to imply dispatching on hypervector type. However, people will want to introduce and compare different operators applied to the same hypervectors, so this issue is orthogonal to hypervector type.

pop, popleft, replace avoid passing the tensor to pop or replace

Add the function hash_table(keys: Tensor, values: Tensor) -> Tensor to the functional module and the class HashTable to the structures module.

• F.hash_table(keys, values)
• S.HashTable(dim)
• S.HashTable().from_tensors(keys, values)
• HashTable[key] -> value
• HashTable[key] = value -> value
• HashTable.get(key) -> value
• HashTable.add(key, value)
• HashTable.set(key, value)
• HashTable.remove(key)
• del HashTable[key]

Where possible, implement the encoding functions in terms of the generic VSA primitives (bind, bundle, etc.) rather than pytorch primitives. This allows for comparing the effect of different specific implementations of the VSA primtives.

See #73 and #74.

We would like to make the package available on Conda under torchhd. I need some help from someone experienced with setting that up because I have not managed so far.

I think it would be nice if we rename the hypervector creation function to remove _hv because I don't think that is necessary since the normal names don't interfere with any functions and the name is still clear.

In addition I think we should alias the hypervector creation and operation functions in the root of the library:

import torchhd

x = torchhd.random(2, 10000)
torchhd.level(2, 10000)
torchhd.circular(2, 10000)

torchhd.bind(x[0], x[1])
torchhd.bundle(x[0], x[1])
torchhd.permute(x[1])

I find it useful to think of VSAs as analog computers for computation on discrete structures. Hypervectors in the VSA correspond to wires in the electronic analog computer. The direction (angle) of the vector corresponds to the label on the wire (i.e. the variable or what it represents) and the magnitude of the vector corresponds to the voltage on the wire (the value of the variable or the degree of support for the represented thing). In this view the direction and magnitude of a vector are conceptually distinct attributes and have different uses.

The usual similarity functions (between two vectors) take the angle and magnitude information and combine it into a single concept of similarity, but I think it is helpful to keep them more separated because they answer different questions.

The angle between two vectors indicates the degree to which they represent the same thing - so this answers questions about what is represented independently of how much support there is for what is being represented. Note that when the angle between two vectors is anything other than quasi-orthogonal this evidence that with very high probability the vectors are in a bundling relationship (because of concentration of measure in hyperdimensional vector spaces). Angular similarity answers questions about the existence of bundling.

Rather than using the actual angle between vectors it is more convenient to use the cosine of the angle between them (because this is equivalent to the dot product of the normalised arguments). Call this function something like similarity.cos.

Long story short - the natural magnitude similarity function is the dot product, which is interpreted as the magnitude of the projection of one argument vector onto the other. This depends on the magnitudes of the two argument vectors and the angle between them. If you have a dictionary of unit magnitude atomic hypervectors and a weighted bundle of them (a multiset), the dot product of each of the dictionary items with the bundle tells you the degree to which that item is present in the bundle (i.e. the scalar weight of that item in the bundle). So dot product is essential to tell you about the degree to which vectors are present. Call this function something like similarity.dot.

In a recurrent VSA circuit a typical approach is to have the solution represented as a weighted bundle (multiset) with the set of items in the bundle fixed over iterations, and the weights associated with each item evolving to indicate the degree of support for that item being in the solution (e.g. see https://www.researchgate.net/profile/Ross-Gayler/publication/310752006_A_distributed_basis_for_analogical_mapping/links/583c01a708aed5c6148cbe23/A-distributed-basis-for-analogical-mapping.pdf).

Hamming distance is isomorphic to dot product similarity with a rescaling of the boolean values to bipolar values.

Based on A Theoretical Perspective on Hyperdimensional Computing we could presumably calculate the ideal threshold so that the user won't have to specify it.

The current implementation of the binary tree requires knowing and providing the path of a node. Look into classical computing APIs of binary trees to see how to design a better interface.

I think we should work on adding a getting started page to the documentation for people unfamiliar with Python or PyTorch since there are also many other researchers that use Matlab or R for instance that could benefit from using Torchhd for their work but we need to make it accessible for them.

functional.cleanup(input: Tensor, memory: Tensor, threshold=0.0) -> Tensor returns the hypervector in memory that is most similar to the input hypervector. If the similarity is less than the threshold an IndexNotFound error should be thrown.

The data structures currently have threshold values which are used to determine whether an element is contained in the sequence for example. We should look into removing support for x in set in favor of set.contains(x) -> float which returns the more informative cosine similarity of x with the set. Similar changes should be made for the other classes.

QPyTorch can be used to simulate low-precision arithmetic maybe we can see if this can be useful for certain applications with the library.

https://readthedocs.org and add functionality to see the documentation of previous versions.

When binding one set with another, it is equivalent to creating a set with all the cross products but more efficient. We should look into adding a function to help with this.

Random selection (randsel) bundling is an alternate implementation of the bundling operator that has the advantages of being conceptually simple and working on every possible hypervector type (see https://rgayler.github.io/VSA_altitude_hold/vsa_basic_operators.html#6_Add_vectors).

It works by elementwise random selection of elements from the operand hypervectors. That is, given operand hypervectors A and B with elements A_i and B_i, each element C_i of the result hypervector C is a copy of A_i or C_i selected at random with equal probability. It delivers a mixture of the operands while preserving the index information.

Viewing bundle as an abstract data type we define bundle in terms of the similarity relationships of the result to the operands - we want the result of bundling to be similar to the operands.
All the usual similarity functions over pairs of hypervectors are effectively averages over the element indices of elementwise similarity. So for the result to be similar to each of the operands it should include some elements from each of the operands.
Because it is just selecting elements at random and not operating on their values it doesn't matter what type the operands have (boolean, integer, real, complex, etc.). The operands and results are all the same type and dimension.

Binary Spatter Code (boolean) hypervectors with majority-rule bundling for two operands is equivalent to randsel bundling. If both operand are the same value (TRUE or FALSE) the result is the same value. This is identical to randomly sampling one of the two operand values. If the operands are different values (TRUE and FALSE) in BSC the tie is broken at random. This is identical to randomly sampling one of the two operand values.

randsel bundling can be generalised to more than two operands (like BSC bundling) and scalar weighting of the operands.
Given k-many operands, just select between the operands operands at random with equal probability.
To do scalar weighted bundling (i.e multiset encoding) there is a finite positive scalar weight associated with each operand, and the elements of each operand are selected with probability proportional to the weight. (Note that the sum of the probabilities across operands sums to one, so if you are thinking of this as a multiset then you need to be aware that the element weights are effectively normalised.)

Note that for more than 2 arguments randsel bundling is not equivalent to BSC bundling because (in the absence of ties) the element result value is the majority element operand value, whereas for randsel bundling (with equal weights) the operand element values appear as result element values with probability reflecting their occurrence across the operands.

Important: The random selection should be random per application of the operator, not a fixed random selection determined at implementation time. A lot of VSA circuits are effectively feed-forward circuits, for these circuits per-operation randomisation is equivalent to at-initialisation randomisation on a single pass through the circuit. However, at-initialisation randomisation does induce correlations between trials, so makes a (slight) difference to statistics calculated over trials. With recurrent VSA circuits the difference between per-operation and at-initialisation is important because the values are passed repeatedly through the same operations. Per-operation randomisation performs the operand mixing without imposing extra structure, whereas at-initialisation randomisation is more like dividing the circuit into multiple parallel fixed circuits without bundling.

(Per-operation randomisation is probably quite expensive in a conventional digital computer and I don't know whether it upsets auto-differentiation in torch, but it's a useful theoretical baseline and might be relatively cheap in unconventional hardware.)

The current auto-formatting workflow fails for PRs from forked repositories. Maybe we should have the formatter workflow only run on the main branch.