Quantization, sketches and estimators

ScalarQuant

Uniform scalar quantization of vector databases at 2, 4 and 8 bits per coordinate, in two families that differ in where the quantization range comes from:

• Global (SQgu2, SQgu4, SQgu8, and GlobalQuantDatabase): a single min/scale pair shared by the whole database, estimated once from a sample of all coordinate values. Codes of different vectors are directly comparable, so the distances consume them with SIMD and nothing is stored per vector.
• Per-vector (SQu2, SQu4, SQu8): every stored vector carries its own min/scale (SQMinC), placed on that vector's own coordinates by a RangePolicy: by default (AutoRange, since 1.6.4) a symmetric mean ± k·σ with a calibrated k at 2 bits, with k searched per vector on a histogram at 4, and the vector's extrema at 8. A vector living on a different scale than the rest is quantized as faithfully as any other; the price is 8 bytes per vector and a distance that must fold both vectors' parameters in (see SQu8's dotu8).

Choosing a family

Use a global quantizer when what you will quantize – the database, the queries, and whatever arrives after the index is built – comes from the same distribution as the sample the range was estimated from. That is the common case for normalized embeddings out of a single model, and it is the cheaper family at query time.

Use a per-vector quantizer when coordinates may move outside what that sample showed: several embedding models or preprocessing pipelines mixed into one database, unnormalized data, scales that drift as the collection grows. A global range is an assumption about the whole database, and a vector that breaks it saturates against the ends of the range, losing precisely the large coordinates distances depend on. A per-vector range cannot be broken, because it is derived from the vector itself.

The trade has a second edge. A per-vector range set on the vector's extrema lets one outlying coordinate stretch it and cost resolution for every other coordinate of that vector, which is why the default policy places a symmetric range at 2 and 4 bits and lets a few coordinates saturate (measured: +0.19-0.20 of recall@10 at 2 bits and +0.006-0.015 at 4 on yahooaq and ccnews; at 8 bits the extrema stay). Per-vector adapts to scale differences between vectors; it does not help against a heavy-tailed coordinate distribution, which costs resolution in both families.

Choosing the range of a global quantizer

The range is a real parameter, and its best value depends on the code width, not only on the data. Clipping at ±L (in units of the coordinate standard deviation) trades two errors against each other: the step is 2L / levels, so a wider range coarsens every code, while a narrower one saturates more mass. More levels make the step cheap, so the optimal L grows with the width. For a Gaussian marginal, minimizing quantization error over sqglobalscale's own reconstruction grid gives

Those factors are specific to this grid, and at 2 bits that matters. sqglobalscale places levels + 1 reconstruction points spanning [min, max] inclusive, 2L / levels apart – not the 2^bits cells of width 2L / 2^bits that textbook tables assume. At 2 bits the two differ by a third in step size (2L/3 against 2L/4), which moves the optimum from the textbook 2.0σ to 1.50σ; at 8 bits they agree to within 0.4% (255 points against 256 cells), which is why the difference is easy to miss.

The default policy here (quant=[0.025, 0.975], i.e. ±1.96σ under a Gaussian) is therefore too wide for SQgu2 and too tight for SQgu8. Measured against the benchmarks' own gold standards on three ANN datasets (600K-3M vectors, 384 dimensions, coordinate kurtosis 3.04-3.05, so Gaussian to measurement precision), moving from the default to the optimal range gained 0.033 to 0.057 of recall@10 at 2 bits and 0.005 to 0.010 at 4 bits; at 8 bits, on five datasets, the same fix was worth 0.027 to 0.079.

Three practical notes on estimating it:

• Prefer passing minmax as an interior quantile times a width-dependent factor over asking for an extreme quantile directly. The default sample holds sqrt(length(X)) values (~15K for a 600K x 384 database), which puts the 8-bit optimum's 99.99th percentile beyond anything the sample can resolve, while the 98th is stable.
• Anchor the range on the median, not on the origin or on the quantile interval's midpoint, and widen each side by its own distance to the median. All three agree on a symmetric centered marginal, and each of the other two breaks on a different shape: scaling both quantiles by one factor displaces a shifted marginal as it widens it (measured at 1.3-1.4x the error of the best range), while scaling about the midpoint drags min below zero on one-signed data such as ReLU outputs, spending levels where there is no data (2.6x the error, -0.046 of recall at 4 bits).
• Never take the range from the extrema. One outlier fixes it for the whole database and the codes collapse onto the middle: over the same datasets, min/max cost 0.15 at 2 bits where the extremes sat at 5.9σ, and 0.67 where a single vector pushed them to 11.7σ. It is also not reproducible, since the sample maximum grows like σ sqrt(2 ln n).
• The factors above assume a near-Gaussian marginal, which a random rotation upstream guarantees and ordinary text embeddings satisfy. Under heavier tails the optimum moves out fast – near 7σ at 8 bits for a Laplace marginal – so when the data has not been Gaussianized, pick L by minimizing quantization error over the sample rather than trusting the table. That search costs milliseconds and tracked the recall optimum closely here, though at 2 bits it lands slightly wide: recall peaked a little tighter than the error did.

SQMinC(min::Float32, c::Float32)

Internal helper struct that stores the per-vector dequantization parameters used by the scalar quantization schemes in ScalarQuant (i.e., SQu2, SQu4, SQu8). Given a quantized (integer) coordinate q, the corresponding approximate original value is recovered as q * c + min.

SQVec{B,VEC<:AbstractVector{UInt8}}

SQVec{B}(v::AbstractVector)
SQVec{B}(E::SQMinC, V::AbstractVector{UInt8})
SQVec{B}(E::SQMinC, V::AbstractVector{UInt8}, Sa, Saa)

A single vector quantized to B bits per coordinate (2, 4 or 8): the packed codes V (codesperbyte(B) coordinates to a byte, low bits first), the affine dequantization parameters E (coordinate i is code * E.c + E.min), and the two code sums Sa = Σ codes and Saa = Σ codes² that let every distance in this module be computed from one integer pass over the codes (see the note above codesums in codes.jl). Indexing (qvec[i]) dequantizes coordinate i to a Float32.

It is the element type of every quantized database here, in both families: a globally quantized vector is a per-vector one whose E happens to be shared with the rest of its database. SQu8Vec, SQu4Vec and SQu2Vec are aliases for SQVec{8}, SQVec{4} and SQVec{2}.

The first constructor quantizes v on its own range (quantvector!, placed by range). length(v) must be a multiple of codesperbyte(B) (2 at 4 bits, 4 at 2 bits), or an ArgumentError is thrown: pad v if needed, and then pad any plain vector later compared against the result the same way, since the mixed distances index it positionally. The other two take stored codes back as they are, recomputing the sums from the codes when they are not given.

quantvector!(::Val{B}, vout::AbstractVector{UInt8}, v::AbstractVector; range=DEFAULT_RANGE, eps=1f-6) -> SQMinC

Quantizes v into vout on its own range, placed by range (a RangePolicy; the default, AutoRange, resolves by width), so a vector on a different scale from the rest is encoded as faithfully as any other. Returns the SQMinC that dequantizes it; eps keeps a constant vector's range from collapsing.

sqdistortion(S, levels, min, max)

Mean squared error the codes of a global quantizer over [min, max] would inflict on the values S, counting both halves of what the range trades: the rounding of what falls inside and the saturation of what falls outside. This is the quantity sqautorange minimizes, and it is computed through the very same arithmetic sqglobalscale hands the quantizers, so the range it picks is optimal for the grid that will actually be used rather than for an idealized one.

The largest code width B can hold, 2^B - 1: 255, 15 or 3.

QuantDatabase{B,P,DB<:AbstractDatabase} <: AbstractDatabase

QuantDatabase{B}(E, Q; dim=nothing, Sa=nothing, Saa=nothing)

An AbstractDatabase of vectors quantized to B bits per coordinate (2, 4 or 8): the packed codes live in any database of UInt8 vectors, Q::DB, next to what dequantizes them, E::P, which is either one SQMinC per vector (P a Vector{SQMinC}: the per-vector family, SQu8Database and its siblings) or a single one shared by the whole database (P == SQMinC: the global family, GlobalQuantDatabase). Indexing yields an SQVec, so every distance in this module applies to both families alike. See the module docstring for how to choose between them, and the aliases above for the constructors that quantize a matrix.

Storage

Q decides what the database can do, exactly as it does for unquantized vectors:

• MatrixDatabase, the default of every constructor that quantizes a matrix: static and fastest;
• BlockMatrixDatabase: grows a block at a time. push_item!/append_items! quantize each incoming vector with the database's own parameters and push its codes, so a SearchGraph can be built over a quantized database one item at a time;
• MMapMatrixDatabase: grows too, on disk, and outlives the process. Reopen the file and hand it back to the constructor with the parameters to get the same database;
• a VectorDatabase of Vector{UInt8}, or anything else indexable: whatever a custom store needs.

A constructor that quantizes a matrix takes storage, a function from the Matrix{UInt8} of codes to the database that will hold them (MatrixDatabase, BlockMatrixDatabase, or a closure that fills an mmap file). The constructor above takes codes back as they are – a matrix, wrapped in a MatrixDatabase, or any database – and an empty database of the right row size is how a growable one starts:

julia> using SimilaritySearch, SimilaritySearch.ScalarQuant

julia> X = randn(Float32, 64, 1000);

julia> db = GlobalQuantDatabase(8, BlockMatrixDatabase(64, UInt8), extrema(X));   # empty

julia> append_items!(db, MatrixDatabase(X)); length(db)
1000

julia> db == GlobalQuantDatabase(8, X; minmax=extrema(X))      # same codes, byte for byte

Growth is serial, like every other database here: concurrent push_item!s need a lock.

Fields

• E, Q: the two that make a database, and the two that have to be persisted;
• dim: the dimension of the vectors, a multiple of the coordinates packed per byte (2 at 4 bits, 4 at 2 bits). A dimension that does not fill its last byte is rejected: pad the vectors, and then pad every plain vector compared against the database the same way, since the mixed distances index it positionally;
• Sa, Saa: per vector, Σ codes and Σ codes², which every distance reads (see the note above codesums in codes.jl). They are sums of the codes, measured when a vector is quantized: sums of the original Float32 vector, though more accurate, lost 0.01 to 0.41 of recall@10 on the SISAP 2025 benchmarks, because quantization rescales each vector slightly and only the code sums carry that rescaling (issue #87). Being derived from the codes they are recomputed on reconstruction unless passed in, which skips one pass over the codes when an mmap-backed database is reopened.

Arguments

• E: a Vector{SQMinC} with one entry per stored vector, or one SQMinC for all
• Q: the codes, as an AbstractMatrix{UInt8} (one column per vector) or an AbstractDatabase of UInt8 vectors

Keyword Arguments

• dim: the vector dimension; read off Q when not given, which an empty VectorDatabase cannot provide
• Sa, Saa: the code sums, when the caller kept them; recomputed from Q otherwise

quantize(db::QuantDatabase, v::AbstractVector) -> SQVec

Quantizes v the way db's vectors are: with db's shared parameters in the global family, so the result is comparable with what it stores; on v's own range (the default policy) in the per-vector one, where db only fixes the dimension. This is what push_item! stores, and what a query must go through to be compared as codes against codes (a Float32 query needs no quantization: the mixed distances take it as it is).

Whether every stored vector shares one min/scale pair (the global family).

The code width of q, in bits.

The code width of the stored vectors, in bits.

SqL2()

The squared Euclidean distance between two quantized vectors (SQVec), or between a quantized vector and a plain vector, avoiding the square root computed by L2.

Between two quantized vectors it never dequantizes: with equal scales it is one exact integer pass over the codes, and otherwise the expansion of Σ (a·cA + mA - b·cB - mB)² into one integer dot product plus the per-vector sums. Against a plain Float32 vector it unpacks the codes to floats with SIMD and accumulates (â - b)².

L2()

The Euclidean ($L_2$) distance between two quantized vectors (SQVec), or between a quantized vector and a plain vector: the square root of SqL2.

L1()

The Manhattan ($L_1$) distance between two quantized vectors (SQVec), or between a quantized vector and a plain one. evaluate dequantizes coordinate by coordinate and accumulates the absolute value of the differences.

NormCosine()

Similar to Dist.NormCosine but for quantized vectors (SQVec): it assumes that the original (pre-quantization) vectors were already normalized, and therefore reduces to one minus the dot product,

\[1 - \sum_i {u_i v_i}\]

Between two quantized vectors the dot product comes from the integer expansion behind Cosine; against a plain vector it is accumulated coordinate by coordinate. When the data may have drifted off the unit sphere through quantization, Cosine corrects for it and ranks better at every width.

Cosine()

Cosine dissimilarity (1 - cos) between two quantized vectors, computed from their codes: the dot product of the dequantized vectors divided by their dequantized norms.

Both halves of that need per-vector quantities the codes alone do not carry, and both are already stored by every SQVec (Sa = Σ codes, Saa = Σ codes²):

• the dot product expands to cA·cB·Σab + cA·mB·Σa + cB·mA·Σb + n·mA·mB, so dropping everything but Σab – as a raw-code dot product does – only preserves order when the offsets are zero. On centered data (any ordinary embedding) it does not: measured over 20k unit-norm vectors in dim 128, ranking by the raw code dot product gave recall@10 of 0.005 against exact cosine, while the full expansion gave 0.97 (issue #77);
• the norms, ‖â‖² = c²·Saa + 2·c·m·Sa + n·m², correct the drift quantization leaves in a vector that was normalized before being quantized. Against a plain vector the query's norm is computed on the spot. That matters more the fewer bits there are: on non-negative data, recall@10 went 0.9685 -> 0.979 at 8 bits, 0.590 -> 0.666 at 4 bits, and 0.111 -> 0.126 at 2 bits.

Both must come from the codes, not from the original vector: on the SISAP 2025 benchmarks, norms measured from the Float32 vector instead lost 0.01-0.05 of recall@10 and sums of the original vector in the dot product lost 0.01-0.41, because quantization rescales each vector slightly and the code sums carry exactly that rescaling (issue #87).

SQu2

Per-vector (per-column) 2-bit scalar quantization: quantize stores four 2-bit codes packed per UInt8, each column keeping its own min/scale computed from its extrema. Accessed as ScalarQuant.SQu2.quantize, etc.

The vector type, the database and the distances are the module-wide SQVec, QuantDatabase and SqL2/L2/L1/ NormCosine; the names here are aliases kept for the per-width API.

quantize(X::AbstractMatrix; storage=MatrixDatabase)

Scalar-quantizes each column (vector) of X to 2 bits per coordinate, four 2-bit codes packed per UInt8. This reduces the memory footprint of a database of vectors by roughly a factor of 16 with respect to Float32 at the cost of precision. Each column is quantized independently using its own minimum and scale factor, computed from the extrema of the column so that its whole range [min, max] maps onto the codes; size(X, 1) must be a multiple of 4.

Returns an SQu2Database, an AbstractDatabase whose db[i] is an SQu2Vec that can be indexed to retrieve dequantized Float32 coordinates; storage chooses the database that holds the codes (see QuantDatabase).

Examples

julia> using SimilaritySearch

julia> X = rand(Float32, 8, 1000);

julia> db = ScalarQuant.SQu2.quantize(X);

julia> db[1][1]  # dequantized approximation of X[1, 1]

quantize(db::SQu2Database, v::AbstractVector)

Quantizes a single vector v to 2 bits per coordinate, the same way as the vectors already stored in db, returning an SQu2Vec. Each vector's own min/scale comes from its own extrema, so this does not read db's parameters; db only fixes the (padded) dimension v must have. It is what push_item! stores, and how a query is quantized to be compared as codes against codes with L1/L2/SqL2/NormCosine.

Examples

julia> using SimilaritySearch

julia> X = rand(Float32, 8, 1000);

julia> db = ScalarQuant.SQu2.quantize(X);

julia> qv = ScalarQuant.SQu2.quantize(db, rand(Float32, 8));  # quantized the same way as db's vectors

SQu2Database(X::AbstractMatrix; storage=MatrixDatabase)
SQu2Database(E::AbstractVector{SQMinC}, Q; dim=nothing, Sa=nothing, Saa=nothing)

The per-vector 2-bit QuantDatabase: four 2-bit codes packed per UInt8, each stored vector carrying its own min/scale pair (E[i]) computed from that vector's own extrema. Indexing yields an SQu2Vec, which the distances consume without dequantizing.

The first constructor quantizes X (it is what quantize calls) into the storage storage builds from the matrix of codes. The second takes the two fields back as they are, quantizing nothing – that is the one to use after reading E/Q from storage, and, with an empty growable Q, the one that starts a database push_item! can grow: see QuantDatabase for the storage choices and the growth interface. The two must agree: exactly one SQMinC per stored vector.

SQu4

Per-vector (per-column) 4-bit scalar quantization: quantize stores two 4-bit codes packed per UInt8, each column keeping its own min/scale computed from its extrema. Accessed as ScalarQuant.SQu4.quantize, etc.

The vector type, the database and the distances are the module-wide SQVec, QuantDatabase and SqL2/L2/L1/ NormCosine; the names here are aliases kept for the per-width API.

quantize(X::AbstractMatrix; storage=MatrixDatabase)

Scalar-quantizes each column (vector) of X to 4 bits per coordinate, two 4-bit codes packed per UInt8. This reduces the memory footprint of a database of vectors by roughly a factor of 8 with respect to Float32 at the cost of precision. Each column is quantized independently using its own minimum and scale factor, computed from the extrema of the column so that its whole range [min, max] maps onto the codes; size(X, 1) must be a multiple of 2.

Returns an SQu4Database, an AbstractDatabase whose db[i] is an SQu4Vec that can be indexed to retrieve dequantized Float32 coordinates; storage chooses the database that holds the codes (see QuantDatabase).

Examples

julia> using SimilaritySearch

julia> X = rand(Float32, 8, 1000);

julia> db = ScalarQuant.SQu4.quantize(X);

julia> db[1][1]  # dequantized approximation of X[1, 1]

quantize(db::SQu4Database, v::AbstractVector)

Quantizes a single vector v to 4 bits per coordinate, the same way as the vectors already stored in db, returning an SQu4Vec. Each vector's own min/scale comes from its own extrema, so this does not read db's parameters; db only fixes the (padded) dimension v must have. It is what push_item! stores, and how a query is quantized to be compared as codes against codes with L1/L2/SqL2/NormCosine.

Examples

julia> using SimilaritySearch

julia> X = rand(Float32, 8, 1000);

julia> db = ScalarQuant.SQu4.quantize(X);

julia> qv = ScalarQuant.SQu4.quantize(db, rand(Float32, 8));  # quantized the same way as db's vectors

SQu4Database(X::AbstractMatrix; storage=MatrixDatabase)
SQu4Database(E::AbstractVector{SQMinC}, Q; dim=nothing, Sa=nothing, Saa=nothing)

The per-vector 4-bit QuantDatabase: two 4-bit codes packed per UInt8, each stored vector carrying its own min/scale pair (E[i]) computed from that vector's own extrema. Indexing yields an SQu4Vec, which the distances consume without dequantizing.

The first constructor quantizes X (it is what quantize calls) into the storage storage builds from the matrix of codes. The second takes the two fields back as they are, quantizing nothing – that is the one to use after reading E/Q from storage, and, with an empty growable Q, the one that starts a database push_item! can grow: see QuantDatabase for the storage choices and the growth interface. The two must agree: exactly one SQMinC per stored vector.

SQu8

Per-vector (per-column) 8-bit scalar quantization: quantize stores one UInt8 code per coordinate, each column keeping its own min/scale computed from its extrema. Accessed as ScalarQuant.SQu8.quantize, etc. See also SQgu8 for a variant that shares a single pair of quantization parameters across all columns.

The vector type, the database and the distances are the module-wide SQVec, QuantDatabase and SqL2/L2/L1/ NormCosine; the names here are aliases kept for the per-width API.

quantize(X::AbstractMatrix; storage=MatrixDatabase)

Scalar-quantizes each column (vector) of X to 8 bits per coordinate, one UInt8 code per coordinate. This reduces the memory footprint of a database of vectors by roughly a factor of 4 with respect to Float32 at the cost of precision. Each column is quantized independently using its own minimum and scale factor, computed from the extrema of the column so that its whole range [min, max] maps onto the codes.

Returns an SQu8Database, an AbstractDatabase whose db[i] is an SQu8Vec that can be indexed to retrieve dequantized Float32 coordinates; storage chooses the database that holds the codes (see QuantDatabase).

Examples

julia> using SimilaritySearch

julia> X = rand(Float32, 8, 1000);

julia> db = ScalarQuant.SQu8.quantize(X);

julia> db[1][1]  # dequantized approximation of X[1, 1]

quantize(db::SQu8Database, v::AbstractVector)

Quantizes a single vector v to 8 bits per coordinate, the same way as the vectors already stored in db, returning an SQu8Vec. Each vector's own min/scale comes from its own extrema, so this does not read db's parameters; db only fixes the (padded) dimension v must have. It is what push_item! stores, and how a query is quantized to be compared as codes against codes with L1/L2/SqL2/NormCosine.

Examples

julia> using SimilaritySearch

julia> X = rand(Float32, 8, 1000);

julia> db = ScalarQuant.SQu8.quantize(X);

julia> qv = ScalarQuant.SQu8.quantize(db, rand(Float32, 8));  # quantized the same way as db's vectors

SQu8Database(X::AbstractMatrix; storage=MatrixDatabase)
SQu8Database(E::AbstractVector{SQMinC}, Q; dim=nothing, Sa=nothing, Saa=nothing)

The per-vector 8-bit QuantDatabase: one UInt8 code per coordinate, each stored vector carrying its own min/scale pair (E[i]) computed from that vector's own extrema. Indexing yields an SQu8Vec, which the distances consume without dequantizing.

The first constructor quantizes X (it is what quantize calls) into the storage storage builds from the matrix of codes. The second takes the two fields back as they are, quantizing nothing – that is the one to use after reading E/Q from storage, and, with an empty growable Q, the one that starts a database push_item! can grow: see QuantDatabase for the storage choices and the growth interface. The two must agree: exactly one SQMinC per stored vector.

2-bit SQVec, four 2-bit codes packed per UInt8. SQu2Vec(v) quantizes v on its own extrema.

4-bit SQVec, two 4-bit codes packed per UInt8. SQu4Vec(v) quantizes v on its own extrema.

8-bit SQVec, one UInt8 code per coordinate. SQu8Vec(v) quantizes v on its own extrema.

sqglobalscale(levels::Integer, min, max)

The scale factor shared by every global quantizer (SQgu2/SQgu4/SQgu8): maps the range [min, max] onto the levels + 1 integer codes 0:levels as code = round(clamp((x - min) * c, 0, levels)). The 1e-6 in the denominator keeps a degenerate (min == max) range from producing Inf.

Arguments

• levels: the largest code the target width can hold (3, 15 or 255)
• min, max: the global value range being mapped

sqautorange(V, levels; samplesize=0, factors=0.30f0:0.05f0:3.0f0, passes=3)

Picks the (min, max) a global quantizer of levels + 1 codes should use for V, by minimizing sqdistortion over a sample of V itself. Returns a (Float32, Float32) pair, ready for minmax.

The search runs over each side separately: the range is anchored at the median and each end is placed at its own multiple of that side's distance to the [0.02, 0.98] quantiles, (min, max) = (med - a(med - q02), med + b(q98 - med)), with a and b found by coordinate descent over factors. Three properties come out of that parametrization, and all three were measured to matter:

• the right width follows the code width. Widening coarsens every code by 2L/levels while narrowing saturates more mass, and more levels make the first cheap, so the optimum moves out as levels grows: for a Gaussian marginal it sits near 1.50σ at 2 bits and 3.93σ at 8. One fixed quantile pair cannot serve both; searching finds each.
• the median anchors it. Scaling both quantiles by a single factor displaces a marginal that is not centered as it widens it, and anchoring at the quantile interval's midpoint instead drags min below zero on one-signed data (ReLU outputs, term weights), spending codes where there is no data. Anchoring at the median does neither, and med == q02 pins a one-signed range at its own floor for free.
• no distributional assumption. Tabulated loading factors are Gaussian values, and the optimum moves out fast under heavier tails (near 7σ at 8 bits for a Laplace marginal). The search reads the shape off the sample instead.

The sample is 1024 (levels + 1) values, or sqrt(length(V)) when that is larger, because the tail the optimum clips shrinks with the code width and has to be visible in the sample to be placed. Cost is a few hundred passes over it – milliseconds even at 8 bits – against a quantization that touches every entry.

sqrange(V, levels; minmax=nothing, quant=nothing, samplesize=0)

The (min, max) every global quantizer in this module resolves before encoding: minmax verbatim when given, else the fixed quantile pair quant when given, else sqautorange's search. Pass quant=[0.025, 0.975] to get the fixed-quantile behaviour that used to be the default.

SQgu2

Global (database-wide) 2-bit scalar quantization: quantize maps every coordinate of every vector using a single shared min/scale pair, packing four 2-bit codes per UInt8, and SqL2 compares the resulting codes directly with SIMD. Accessed as ScalarQuant.SQgu2.quantize, etc.

It is the coarsest member of the global family (SQgu2/SQgu4/SQgu8): four codes fit in one byte, so a Float32 database shrinks by 16x, at the cost of only four levels per coordinate. Codes are laid out low bits first (coordinate 4k+1 in bits 0:1, 4k+2 in bits 2:3, and so on), exactly like SQu2's per-column variant.

quantize(X::AbstractMatrix; minmax=nothing, quant=nothing, samplesize=0)

Scalar-quantizes every entry of X to 2 bits using a single, global pair of dequantization parameters shared by all columns, unlike SQu2's quantize which computes an independent min/scale per column. As with SQgu4/SQgu8, this is useful when the columns of X share a comparable value range, since a single global range provides enough precision while being cheaper to compute and store.

Codes are packed four per UInt8 (low bits first), so the returned matrix has cld(size(X, 1), 4) rows. Packing four dimensions into a single byte, combined with a global (rather than per-column) min/scale, lets SqL2 and SqL2 operates directly on the packed codes with SIMD, without any per-element dequantization: since every column shares the same affine mapping, comparisons and (squared) differences computed in code space are already proportional to the ones in the original space.

The global [min, max] range is estimated by sampling entries of X and taking quantiles of the sample (to be robust to outliers), unless it is provided explicitly via minmax. Every entry x is then mapped as round(clamp((x - min) * c, 0, 3)) with c = 3 / (max - min + 1e-6).

Arguments

• X: the matrix to quantize; each entry is quantized independently but using shared min/max values
• minmax: an optional (min, max) tuple giving the value range to use; when nothing (the default) the range is chosen from a random sample of the entries of X by sqautorange, or by quant when that is given
• quant: a fixed lower/upper quantile pair (of the sampled entries of X) to use as the range instead of searching for it. nothing (the default) runs sqautorange, which places each end of the range where it minimizes the error the codes would incur – the optimum moves with the code width, so a fixed pair cannot be right at 2, 4 and 8 bits at once. Pass [0.025, 0.975] for the pre-search behaviour
• samplesize: the number of entries sampled (with replacement) from X to estimate the quantiles; when 0 (the default) it is set to ceil(Int, length(X)^0.5)

Examples

julia> using SimilaritySearch

julia> X = rand(Float32, 8, 1000);

julia> Q = ScalarQuant.SQgu2.quantize(X; minmax=(0f0, 1f0));  # explicit range

julia> size(Q), eltype(Q)  # (2, 1000), UInt8

quantize(v::AbstractVector; minmax=nothing, quant=nothing, samplesize=0)

Scalar-quantizes a single vector v to 2 bits, using the same global scheme as quantize(X::AbstractMatrix), producing a Vector{UInt8} (four codes per byte) of length cld(length(v), 4), instead of a Matrix{UInt8}.

Arguments

• v: the vector to quantize
• minmax: an optional (min, max) tuple giving the value range to use; when nothing (the default) the range is estimated from a random sample of v's entries using quant. Must match the dataset's minmax if v is to be compared against an existing quantized dataset.
• quant: a fixed lower/upper quantile pair to use as the range instead of searching for it; nothing (the default) runs sqautorange. Used when minmax is not given
• samplesize: the number of entries sampled (with replacement) from v to estimate the quantiles; when 0 (the default) it is set to ceil(Int, length(v)^0.5)

quantize!(vout::AbstractVector{UInt8}, v::AbstractVector, minmax) -> vout

In-place, allocation-free quantize of a single vector into a caller-provided vout of length cld(length(v), 4), using the explicit (min, max) range minmax. Intended for encoding loops that reuse their output buffer (see Projections.quantsketch); the range is never estimated here, precisely so every vector encoded through it stays comparable.

SqL2()

Squared Euclidean distance between two vectors quantized with quantize (four globally-scaled 2-bit codes per byte). Since both vectors share the same global min/scale, the squared difference of the raw codes is proportional to the squared difference of the original values, so evaluate accumulates squared code differences directly, unpacking each byte's four 2-bit fields with SIMD, without any per-element dequantization.

SQgu4

Global (database-wide) 4-bit scalar quantization: quantize maps every coordinate of every vector using a single shared min/scale pair, packing two 4-bit codes per UInt8, and SqL2 compares the resulting codes directly with SIMD. Accessed as ScalarQuant.SQgu4.quantize, etc.

quantize(X::AbstractMatrix; minmax=nothing, quant=nothing, samplesize=0)

Scalar-quantizes every entry of X to 4 bits using a single, global pair of dequantization parameters shared by all columns, unlike SQu4's quantize which computes an independent min/scale per column. As with SQgu8's quantize, this is useful when the columns of X share a comparable value range, since a single global range provides enough precision while being cheaper to compute and store.

Codes are packed two per UInt8 (low nibble, high nibble), exactly like SQu4's, so the returned matrix has ceil(Int, size(X, 1) / 2) rows. Packing pairs of dimensions into a single byte, combined with a global (rather than per-column) min/scale, lets SqL2 operates directly on the packed codes with SIMD, without any per-element dequantization: since every column shares the same affine mapping, comparisons and (squared) differences computed in code space are already proportional to the ones in the original space.

The global [min, max] range is estimated by sampling entries of X and taking quantiles of the sample (to be robust to outliers), unless it is provided explicitly via minmax. Every entry x is then mapped as round(clamp((x - min) * c, 0, 15)) with c = 15 / (max - min + 1e-6).

Arguments

• X: the matrix to quantize; each entry is quantized independently but using shared min/max values
• minmax: an optional (min, max) tuple giving the value range to use; when nothing (the default) the range is chosen from a random sample of the entries of X by sqautorange, or by quant when that is given
• quant: a fixed lower/upper quantile pair (of the sampled entries of X) to use as the range instead of searching for it. nothing (the default) runs sqautorange, which places each end of the range where it minimizes the error the codes would incur – the optimum moves with the code width, so a fixed pair cannot be right at 2, 4 and 8 bits at once. Pass [0.025, 0.975] for the pre-search behaviour
• samplesize: the number of entries sampled (with replacement) from X to estimate the quantiles; when 0 (the default) it is set to ceil(Int, length(X)^0.5)

Examples

julia> using SimilaritySearch

julia> X = rand(Float32, 8, 1000);

julia> Q = ScalarQuant.SQgu4.quantize(X; minmax=(0f0, 1f0));  # explicit range

julia> size(Q), eltype(Q)  # (4, 1000), UInt8

quantize(v::AbstractVector; minmax=nothing, quant=nothing, samplesize=0)

Scalar-quantizes a single vector v to 4 bits, using the same global scheme as quantize(X::AbstractMatrix), producing a Vector{UInt8} (nibble-packed, two codes per byte) of length ceil(Int, length(v) / 2), instead of a Matrix{UInt8}.

Arguments

• v: the vector to quantize
• minmax: an optional (min, max) tuple giving the value range to use; when nothing (the default) the range is estimated from a random sample of v's entries using quant. Must match the dataset's minmax if v is to be compared against an existing quantized dataset.
• quant: a fixed lower/upper quantile pair (of the sampled entries of v) to use as the range instead of searching for it. nothing (the default) runs sqautorange, which places each end of the range where it minimizes the error the codes would incur – the optimum moves with the code width, so a fixed pair cannot be right at 2, 4 and 8 bits at once. Pass [0.025, 0.975] for the pre-search behaviour
• samplesize: the number of entries sampled (with replacement) from v to estimate the quantiles; when 0 (the default) it is set to ceil(Int, length(v)^0.5)

Examples

julia> using SimilaritySearch

julia> minmax = (0f0, 1f0);

julia> X = rand(Float32, 8, 1000);

julia> Q = ScalarQuant.SQgu4.quantize(X; minmax);  # dataset, using an explicit range

julia> q = rand(Float32, 8);

julia> qv = ScalarQuant.SQgu4.quantize(q; minmax);  # query, using the *same* range

julia> length(qv), eltype(qv)  # (4, UInt8)

quantize!(vout::AbstractVector{UInt8}, v::AbstractVector, minmax) -> vout

In-place, allocation-free quantize of a single vector into a caller-provided vout of length cld(length(v), 2), using the explicit (min, max) range minmax. Intended for encoding loops that reuse their output buffer (see Projections.quantsketch); the range is never estimated here, precisely so every vector encoded through it stays comparable.

SqL2()

Squared Euclidean distance between two vectors quantized with quantize (nibble-packed, globally-scaled 4-bit codes). Since both vectors share the same global min/scale, the squared difference of the raw codes is proportional to the squared difference of the original values, so evaluate accumulates squared code differences directly, unpacking each byte's low and high nibble with SIMD, without any per-element dequantization.

SQgu8

Global (database-wide) 8-bit scalar quantization: quantize maps every coordinate of every vector using a single shared min/scale pair, and SqL2 compares the resulting codes directly with SIMD. Accessed as ScalarQuant.SQgu8.quantize, etc.

quantize(X::AbstractMatrix; minmax=nothing, quant=nothing, samplesize=0)

Scalar-quantizes every entry of X to 8 bits (UInt8) using a single, global pair of dequantization parameters shared by all columns, unlike SQu8's quantize which computes an independent min/scale per column. This is useful, e.g., when the columns of X are known to share a comparable value range and a single global range provides enough precision while being cheaper to compute and store.

The global [min, max] range is estimated by sampling entries of X and taking quantiles of the sample (to be robust to outliers), unless it is provided explicitly via minmax. Every entry x is then mapped as round(clamp((x - min) * c, 0, 255)) with c = 255 / (max - min + 1e-6).

Arguments

• X: the matrix to quantize; each entry is quantized independently but using shared min/max values
• minmax: an optional (min, max) tuple giving the value range to use; when nothing (the default) the range is chosen from a random sample of the entries of X by sqautorange, or by quant when that is given
• quant: a fixed lower/upper quantile pair (of the sampled entries of X) to use as the range instead of searching for it. nothing (the default) runs sqautorange, which places each end of the range where it minimizes the error the codes would incur – the optimum moves with the code width, so a fixed pair cannot be right at 2, 4 and 8 bits at once. Pass [0.025, 0.975] for the pre-search behaviour
• samplesize: the number of entries sampled (with replacement) from X to estimate the quantiles; when 0 (the default) it is set to ceil(Int, length(X)^0.5)

Examples

julia> using SimilaritySearch

julia> X = rand(Float32, 8, 1000);

julia> Q = ScalarQuant.SQgu8.quantize(X; minmax=(0f0, 1f0));  # explicit range

julia> size(Q), eltype(Q)  # (8, 1000), UInt8

quantize(v::AbstractVector; minmax=nothing, quant=nothing, samplesize=0)

Scalar-quantizes a single vector v to 8 bits (UInt8), using the same global scheme as quantize(X::AbstractMatrix), producing a Vector{UInt8} (one code per coordinate) instead of a Matrix{UInt8}.

Arguments

• v: the vector to quantize
• minmax: an optional (min, max) tuple giving the value range to use; when nothing (the default) the range is estimated from a random sample of v's entries using quant. Must match the dataset's minmax if v is to be compared against an existing quantized dataset.
• quant: a fixed lower/upper quantile pair (of the sampled entries of v) to use as the range instead of searching for it. nothing (the default) runs sqautorange, which places each end of the range where it minimizes the error the codes would incur – the optimum moves with the code width, so a fixed pair cannot be right at 2, 4 and 8 bits at once. Pass [0.025, 0.975] for the pre-search behaviour
• samplesize: the number of entries sampled (with replacement) from v to estimate the quantiles; when 0 (the default) it is set to ceil(Int, length(v)^0.5)

Examples

julia> using SimilaritySearch

julia> minmax = (0f0, 1f0);

julia> X = rand(Float32, 8, 1000);

julia> Q = ScalarQuant.SQgu8.quantize(X; minmax);  # dataset, using an explicit range

julia> q = rand(Float32, 8);

julia> qv = ScalarQuant.SQgu8.quantize(q; minmax);  # query, using the *same* range

julia> length(qv), eltype(qv)  # (8, UInt8)

quantize!(vout::AbstractVector{UInt8}, v::AbstractVector, minmax) -> vout

In-place, allocation-free quantize of a single vector into a caller-provided vout of length length(v), using the explicit (min, max) range minmax. Intended for encoding loops that reuse their output buffer (see Projections.quantsketch); the range is never estimated here, precisely so every vector encoded through it stays comparable.

SqL2()

Squared Euclidean distance between two vectors quantized with quantize (globally-scaled 8-bit codes). Since both vectors share the same global min/scale, the squared difference of the raw codes is proportional to the squared difference of the original values, so evaluate accumulates squared code differences directly with SIMD, widening each UInt8 code to Int32 to safely represent negative differences, without any per-element dequantization.

GlobalQuantDatabase{B}

GlobalQuantDatabase(bits::Integer, X::AbstractMatrix; minmax=nothing, storage=MatrixDatabase, quant=nothing, samplesize=0)
GlobalQuantDatabase(bits::Integer, Q, minmax; dim=nothing, Sa=nothing, Saa=nothing)

The QuantDatabase of vectors quantized to bits bits per coordinate (2, 4 or 8) under one min/scale pair shared by the whole dataset, keeping that pair – and the per-vector code sums – alongside the codes.

That is the difference from calling SQgu8.quantize directly, which hands back a bare Matrix{UInt8} and leaves minmax to the caller. Without the parameters a stored matrix of codes cannot be dequantized at all, so it can only ever be compared against other codes from the same run; with them, a query may stay in its original Float32 form, and a vector arriving later is quantized with the very same range (push_item!), which is what keeps its codes comparable.

Indexing yields an SQVec of the same width the per-vector quantizers produce – a globally quantized vector is a per-vector one whose scale happens to be shared – so every distance in this module works here unchanged, and each takes the path it should:

• SqL2/L2 between two stored vectors hit the equal-scale branch, which is an exact integer pass over the codes;
• SqL2/L2/L1/NormCosine against a plain Float32 vector take the mixed kernels;
• Cosine uses the stored sums for both the dot product's offset terms and the norms (issue #77).

Arguments

• bits: 2, 4 or 8
• X: the matrix to quantize, one column per vector; size(X, 1) must be a multiple of the coordinates packed per byte (2 at 4 bits, 4 at 2 bits), or an ArgumentError is thrown – pad X, and then pad every plain vector compared against the database the same way
• Q, minmax: already-quantized codes (a matrix, or any AbstractDatabase of UInt8 vectors, an empty one to start a growable database) and the pair they were produced with

Keyword Arguments

• minmax: the (min, max) pair to quantize with; estimated from a sample of X by sqrange when not given, exactly as the underlying quantize does
• storage: the database the codes go to, as a function of their Matrix{UInt8}
• dim, Sa, Saa: see QuantDatabase

SQEncoder(quant::Module, X::AbstractMatrix; dist=ScalarQuant.SqL2(), samplesize=4096, rng, range=AutoRange())
SQEncoder(quant::Module, dim::Integer; dist=ScalarQuant.SqL2(), minmax=nothing, range=AutoRange())
SQEncoder(quant::Module, rotation, X::AbstractMatrix; ...)
SQEncoder(quant::Module, rotation, dim::Integer; ...)

The scalar quantizers as the encoder of an AsymmetricSearchGraph: an object is quantized once (after an optional rotation) and stored as its codes, a raw query is kept in Float32 (rotated the same way), and this module's distances evaluate one against the other. It goes through the AbstractEstimator interface – encode, encodequery, evaluate – because that is what the graph navigates with, but it is a codification, not an estimate with an error to exploit: nothing is bounded and nothing is ever re-evaluated.

quant is one of ScalarQuant's quantizer modules, which names the family and the width at once: SQgu8, SQgu4, SQgu2 (one global range for every code, estimated by sqautorange on a sample of X's rotated coordinates, or given as minmax) or SQu8, SQu4, SQu2 (each code with its own range, placed by the range policy – AutoRange by default, resolved by width and calibrated on the sample when one is given, see RangePolicy).

rotation is the object that rotates, a SimilaritySearch.Projections.Rotation – Projections.qr(dim, dim) or Projections.RandomizedHadamard(dim) – or nothing, which quantizes the coordinates as they are and is the default: the forms without a rotation argument rotate nothing. A rotation gives every coordinate the same scale, which is what a single global range rests on; on data whose coordinates already share one (normalized embeddings: on SISAP 2025 ccnews a QR rotation moved recall@10 by less than 0.01 at every width) it changes nothing and costs its flops per query and per inserted item. Rotate when the coordinates' scales are uneven and the global family is wanted anyway; the per-vector family is the other remedy for uneven scales.

encode rotates an object once and quantizes the rotated vector, encodequery rotates the query once and prepares it (an SQQuery: the rotated Float32 vector with its sums and an integer image on its own range), and evaluate is dist, one of this module's (SqL2, L2, L1, NormCosine or Cosine): its factored kernel for the prepared query against a stored code, one integer dot product per pair, and its integer kernels between two stored codes. Nothing is ever rotated or prepared inside an evaluation. sqcodes(e) is the storage to build an AsymmetricSearchGraph over: a QuantDatabase with the estimator's own parameters, whose codes live in dense blocks, and which takes the SQVecs encode produces without re-quantizing them.

RangePolicy

Where the per-vector family places a vector's own range [min, max] before mapping it onto the codes 0:levels(B). AutoRange (the default) picks by width among ExtremaRange, FixedRange, CalibratedRange, HistogramRange, RefinedRange and ExactRange. A policy is applied with vectorrange(policy, v, L) -> (min, step); one that learns from data (the calibrated k) is fitted with calibrate(policy, X, L) when an SQEncoder is built from a matrix, and the encoder keeps the fitted policy. A per-vector database built from a matrix only resolves the policy: it must hold the same codes that growing it with push_item! would.

AutoRange()

The per-vector family's default, resolved by width (resolverange): CalibratedRange at 2 bits, HistogramRange at 4, ExtremaRange at 8.

ExtremaRange()

The vector's own extrema: nothing saturates, and one outlying coordinate stretches the range and coarsens every other code of the vector. The rule at every width up to 1.6.3, and AutoRange's at 8 bits, where a clipped coordinate costs more than the finer step buys back.

FixedRange(k)

mean ± k·σ with a fixed k. The best fixed factors on near-Gaussian coordinates are about 1.5, 2.55 and 4.0 at 2, 4 and 8 bits (defaultk); at 4 bits a fixed factor gives up what a per-vector search gains.

CalibratedRange(; k=0, samplesize=4096, grid=0.5:0.05:4.5)

mean ± k·σ with one k for every vector, calibrated on samplesize vectors of the data by the mean of their own distortions over grid (calibrate, ~0.4 s at 384 dimensions); per vector only the mean and σ are computed, one pass. Uncalibrated (k == 0, e.g. a lone SQVec{B}(v)), it uses defaultk. AutoRange's policy at 2 bits: the same recall as a per-vector histogram search at a third of its cost.

HistogramRange(; bins=64, kmax=4.5)

mean ± k·σ with k searched per vector on a histogram of the deviations |x - mean|/σ, bins bins up to kmax·σ, filled in one pass; each candidate (a bin edge from 0.5) is scored in constant time with the saturation error summed exactly from the histogram and the rounding error modelled as step²/12. About 3 µs per vector at 384 dimensions. AutoRange's policy at 4 bits. 32 and 16 bins were measured below 64 at every width for 0.1-0.2 µs less.

RefinedRange(inner=CalibratedRange(); step=0.1)

The symmetric policy inner's k, then the real grid's sqdistortion of the vector at k - step, k and k + step, keeping the best: three passes more. Refining the calibrated k (4.4 µs) or the histogram's (6.3 µs) recovers most of the exact search at 2 bits.

ExactRange(; kmin=0.5, kmax=4.5)

mean ± k·σ with k minimizing the vector's own sqdistortion: a sweep every 0.25, then every 0.05 around its best (28 passes over the vector, ~32 µs). The reference the others are measured against.

resolverange(policy, B) -> RangePolicy

The concrete policy policy stands for at B bits: AutoRange() becomes the policy chosen for that width, any other policy itself.

calibrate(policy, X::AbstractMatrix, L) -> RangePolicy

The policy fitted to the columns of X at L levels: a CalibratedRange without k gets the one minimizing the mean distortion of samplesize columns; a RefinedRange calibrates its inner policy; every other policy is returned as it is.

defaultk(L) -> Float32

The factor a calibrated policy uses before calibration: the calibrated k measured on yahooaq and ccnews (1.5, 2.55, 4.0 at 3, 15, 255 levels), which are also the Gaussian optima of this grid.

SQQuery{B}(q::AbstractVector) <: AbstractVector{Float32}

A query prepared once for the mixed kernels against SQVec{B} codes (issue #110): what SQEncoder's encodequery returns, and what an AsymmetricSearchGraph evaluates its distance with, once per query and once per inserted item.

With x̂ᵢ = c·aᵢ + m the squared Euclidean distance expands as

Σ x̂ᵢ² − 2 Σ x̂ᵢ qᵢ + ‖q‖²  =  (c²·Saa + 2cm·Sa + n·m²)  −  2·(c·Σ aᵢqᵢ + m·Σqᵢ)  +  ‖q‖²

where the first term comes from the sums every SQVec stores, the last two from this object, and only Σ aᵢqᵢ is evaluated per pair. That dot product is taken with an integer image of the query: qᵢ ≈ mq + sq·uᵢ on the query's own range, uᵢ at 15 bits against 8-bit codes and at 8 bits against 4- and 2-bit codes, stored signed (dᵢ = uᵢ − half) in the planes the packed codes unpack into, so Σ aᵢqᵢ = mq·Sa + sq·(Σ aᵢdᵢ + half·Sa). Measured on SISAP 2025 ccnews against the exact distance: the 8-bit image against 4-bit codes is as accurate as the Float32 query (mean error 0.0165 for both, against 0.0268 for a query quantized on the codes' global range), and the 15-bit image against 8-bit codes is within 1e-4 of it. The kernels are the integer dot products of codes.jl, the symmetric kernels' cost.

It is still an AbstractVector{Float32} – the rotated query, coordinate by coordinate – so everything that takes a plain query keeps working on it: L1, the generic scans, a QuantDatabase quantizing it for the symmetric graph. length(q) must be a multiple of the codes per byte (2 at 4 bits, 4 at 2 bits), as for SQVec.

sqcodes(e::SQEncoder) -> QuantDatabase

An empty, growable QuantDatabase with e's own parameters and dense code blocks: the storage of an AsymmetricSearchGraph over e. It takes the SQVecs encode produces as they are.

sqcodes(e::SQEncoder, X::AbstractMatrix; minbatch=4) -> QuantDatabase

Encodes every column of X in parallel into the storage above.

quantizer(x) -> Module

The quantizer module x was built with, one of SQgu8, SQgu4, SQgu2, SQu8, SQu4, SQu2: the family and the width in one name. Defined for a QuantDatabase, a SQEncoder and RaBitQ.RaBitQVectorFallback.

RaBitQ

The RaBitQ estimator (Gao & Long, 2024) as an AbstractEstimator for an AsymmetricSearchGraph: every object is stored as the sign bits of its randomly rotated vector plus three scalars, and a raw query is evaluated against those bits with an unbiased estimate of the inner product whose error bound is known per object.

est = RaBitQCosine(Projections.qr(dim, dim))    # or RaBitQL2; RandomizedHadamard(dim) rotates in dim log dim
G = AsymmetricSearchGraph(est, rabitqcodes(est))  # a growable storage of codes
append_items!(G, ctx, MatrixDatabase(X))        # raw vectors in, codes stored
search(G, ctx, q, knnqueue(KnnSorted, 10))      # raw query in, rotated once inside

The graph rotates each query once (encodequery), encodes each inserted item once (encode), and evaluates the estimator between the rotated query and the stored bits; it also evaluates two stored codes against each other, which its neighborhood filters need, with the SimHash estimate over the Hamming distance of the bits.

RaBitQCosine(rotation)
RaBitQL2(rotation)

The RaBitQ estimators of the cosine dissimilarity (1 - cos) and of the Euclidean distance between a raw query and a stored code, over rotation, a SimilaritySearch.Projections.Rotation: Projections.qr(dim, dim) or RandomizedHadamard(dim). The dimension is read off it. The rotation is not optional here: the estimate is unbiased and its error bound holds because the sign vector is taken in a uniformly random basis.

Both are AbstractEstimators: encode(est, o) gives the RaBitQCode the storage receives, encodequery(est, q) the RaBitQQuery a raw query becomes, and evaluate(est, q::RaBitQQuery, o::RaBitQCode) the estimated distance. Between two codes, evaluate(est, a::RaBitQCode, b::RaBitQCode) estimates the cosine from the Hamming distance of the bits, cos(π h / dim), which is what the graph's neighborhood filters use.

RaBitQCode{V<:AbstractVector{UInt64}}

What is stored per object: the sign bits of its rotated vector (dim bits in cld(dim, 64) words, low bit first) and three scalars the estimator needs – c = ⟨ō, o⟩ / ‖o‖, the projection of the unit vector onto its own sign vector, which normalizes the estimate; norm = ‖o‖, for the Euclidean distance; and err, the scale of the estimate's confidence interval, sqrt((1 - c²) / c²) · 1.9 / sqrt(dim - 1).

RaBitQQuery{V<:AbstractVector{Float32}}

A raw query prepared once for the estimator: rotated and normalized to a unit vector (r), with the raw norm (norm) kept for the Euclidean distance.

rabitqcodes(est) -> VectorDatabase

An empty, growable storage of the codes est produces, for an AsymmetricSearchGraph over it: RaBitQCodes for the bit estimators, (RaBitQCode, fine) tuples for a RaBitQRefined.

rabitqcodes(e::AbstractRaBitQ, X::AbstractMatrix; minbatch=4) -> VectorDatabase

Encodes every column of X in parallel into a growable storage of codes.

estimatecos(e::AbstractRaBitQ, q::RaBitQQuery, o::RaBitQCode) -> Float32

The RaBitQ estimate of the cosine between the raw query behind q and the object behind o: ⟨ō, q⟩ / (⟨ō, o⟩ ‖q‖), unbiased, with a confidence interval of half-width about errorbound(e, o) (for unit vectors).

errorbound(e, o::RaBitQCode) -> Float32

The half-width of the estimate's confidence interval for the object behind o, sqrt((1 - c²) / c²) · ε₀ / sqrt(dim - 1) with ε₀ = 1.9: the true cosine lies within it of the estimate with probability about 1 - exp(-ε₀² / 2) ≈ 0.84 (measured at 0.80 on unit Gaussian vectors in 384-d). Larger for objects their sign vector represents worse.

RaBitQRefined(coarse::RaBitQCosine, fine::AbstractFallback; τ=Inf32)

Two levels in one estimator: the bits of coarse navigate, and the finer representation beside them re-evaluates, inside the same evaluate, every object whose estimated cosine distance minus its error bound is at or below τ – the objects the bits cannot rule out as far. τ = Inf re-evaluates everything (the fine level decides every distance); a finite τ on the scale of the k-th neighbor's distance, which refinethreshold reads off a sample, keeps the fine evaluations to the candidates that can matter. Between two stored objects (the graph's neighborhood filters) only the bits are compared.

The stored object is a tuple (RaBitQCode, fine); rabitqcodes(e) gives an empty growable storage of the right element type.

abstract type AbstractFallback end

The finer representation a RaBitQRefined keeps beside the bits, in the rotated space: finecode(f, r) builds it from the rotated vector r, and finecos(f, fine, qr) is the cosine between it and the rotated unit query qr. Two are provided: RaBitQExactFallback and RaBitQVectorFallback.

RaBitQExactFallback{T}()   (T = Float32 by default, or Float16)

The rotated vector itself, normalized, in T: 4·dim (2·dim) bytes per object and the exact cosine when it is consulted. This is what RaBitQ's own search re-ranks with, and the case where the bits pay: navigation costs the bit estimate, and the exact evaluation is only paid for the candidates the error bound cannot rule out.

RaBitQVectorFallback(quant::Module, est::AbstractRaBitQ, X::AbstractMatrix; samplesize=4096, rng)
RaBitQVectorFallback(quant::Module, est::AbstractRaBitQ; minmax=nothing)

The rotated unit vector scalar-quantized with quant, one of ScalarQuant's quantizer modules, which names the family and the width at once: SQgu8, SQgu4, SQgu2 (one global range for every code, estimated by sqautorange on a sample of X's rotated unit vectors, or given as minmax) or SQu8, SQu4, SQu2 (each code with its own range, placed by range, a ScalarQuant.RangePolicy; AutoRange() by default, calibrated on the sample when given). est supplies the rotation and the dimension. dim·B/8 + 8 bytes per object, consulted through ScalarQuant's mixed kernels with the dequantized norm.

refinethreshold(est::AbstractRaBitQ, X, k; samplesize=2000, numqueries=200, q=0.9, rng) -> Float32

A τ for RaBitQRefined on the scale of the data: the q-quantile, over numqueries objects of X, of the exact cosine distance to their k-th nearest neighbor within a sample of samplesize objects. A candidate whose estimated distance minus its error bound is within that is one the search may keep, and is worth the fine evaluation.

RandomProjections(map::M) where {M<:AbstractMatrix}

Wraps a projection matrix map of size (indim, outdim) used to reduce the dimension of vectors from indim to outdim via a linear projection (v -> map' * v). This is a standard dimensionality-reduction technique for similarity search: projecting onto a lower-dimensional space reduces both memory usage and distance-computation cost, while approximately preserving relative distances (see the Johnson-Lindenstrauss lemma).

Use gaussian or qr to build a RandomProjections map, and transform/transform! to apply it to vectors or matrices of vectors.

Arguments

• map: the projection matrix, with indim rows and outdim columns

Examples

julia> using SimilaritySearch

julia> rp = SimilaritySearch.Projections.gaussian(128, 32);  # random gaussian projection 128 -> 32

julia> rp2 = SimilaritySearch.Projections.qr(128, 32);  # QR-orthogonalized projection 128 -> 32

gaussian(rng::AbstractRNG, FloatType::Type, indim::Int, outdim::Int)
gaussian(indim::Int, outdim::Int=indim)

Builds a RandomProjections whose map is a dense indim × outdim matrix with entries drawn independently from a Normal distribution with mean zero and standard deviation 1/outdim, whose columns are then normalized to unit norm. This is a Gaussian random projection: unlike qr, the resulting columns are not orthogonal to each other, but generating and applying it is cheaper.

Arguments

• rng: the random number generator to use (defaults to Random.default_rng())
• FloatType: the floating point type of the projection matrix (defaults to Float32)
• indim: the dimension of the input vectors
• outdim: the dimension of the projected vectors (defaults to indim)

Examples

julia> using SimilaritySearch

julia> rp = SimilaritySearch.Projections.gaussian(128, 32);

julia> size(SimilaritySearch.Projections.getmap(rp))
(128, 32)

qr(rng::AbstractRNG, FloatType::Type, indim::Int, outdim::Int)
qr(indim::Int, outdim::Int=indim)

Builds a RandomProjections whose map is the (first outdim columns of the) Q factor of the QR decomposition of a random indim × indim matrix. Unlike gaussian, the resulting projection matrix has orthonormal columns, which makes the projection an isometry up to the subspace it projects onto (distances between projected vectors are not shrunk by non-orthogonality), at the extra cost of computing the QR factorization.

Arguments

• rng: the random number generator to use (defaults to Random.default_rng())
• FloatType: the floating point type of the projection matrix (defaults to Float32)
• indim: the dimension of the input vectors
• outdim: the dimension of the projected vectors (defaults to indim)

outdim(rp::RandomProjections)

Returns the output dimension of the projection rp, i.e., the dimension of the vectors produced by transform/transform!.

outdim(hp::HadamardProjection)

Returns the output dimension of the projection hp, i.e., the dimension of the vectors produced by transform/transform!. Always equal to indim(hp), since HadamardProjection does not reduce dimensionality.

outdim(p::PCAProjection)

Returns the output dimension of the projection p, i.e., the dimension of the vectors produced by transform/transform!.

outdim(m::DistantHyperplanes)

Returns the number of output bits (hyperplanes) of m.

outdim(m::AnchoredDistantHyperplanes)

Returns the number of output bits (hyperplanes) of m.

outdim(B::RandomHyperplanes)

Returns the number of output bits (reference pairs) of B.

outdim(qs::QuantSketch)

The number of components (hyperplanes/directions) of the underlying model, i.e. how many codes each sketch holds – not how many machine words it occupies (see sketchsize).

indim(rp::RandomProjections)

Returns the input dimension of the projection rp, i.e., the dimension that vectors passed to transform/transform! are expected to have.

indim(hp::HadamardProjection)

Returns the input dimension of the projection hp, i.e., the dimension that vectors passed to transform/transform! are expected to have.

indim(p::PCAProjection)

Returns the input dimension of the projection p, i.e., the dimension that vectors passed to transform/transform! are expected to have.

transform(rp::RandomProjections, v::AbstractVector)

Projects the vector v (of length indim(rp)) using rp, returning a new vector of length outdim(rp). Each output coordinate is the dot product of v with the corresponding column of the projection map.

Arguments

• rp: the projection to apply
• v: the input vector to project

transform(rp::RandomProjections, X::AbstractMatrix; minbatch::Int=4)

Projects every column (vector) of X using rp, returning a new matrix with outdim(rp) rows and the same number of columns as X. Columns are projected in parallel using @BATCHES.

Arguments

• rp: the projection to apply
• X: a matrix whose columns are the vectors to project, each of length indim(rp)
• minbatch: minimum number of columns processed per parallel task (see @BATCHES)

Examples

julia> using SimilaritySearch

julia> X = rand(Float32, 128, 1000);

julia> rp = SimilaritySearch.Projections.gaussian(128, 32);

julia> Y = SimilaritySearch.Projections.transform(rp, X);

julia> size(Y)
(32, 1000)

transform(hp::HadamardProjection, v::AbstractVector)

Projects the vector v (of length indim(hp)) using hp, returning a new vector of the same length. Computed as the fast Walsh-Hadamard transform of v (sequency-ordered).

Arguments

• hp: the projection to apply
• v: the input vector to project

transform(hp::HadamardProjection, X::AbstractMatrix; minbatch::Int=4)

Projects every column (vector) of X using hp, returning a new matrix of the same size as X: one fwht! per column, in parallel (see transform!).

Arguments

• hp: the projection to apply
• X: a matrix whose columns are the vectors to project, each of length indim(hp)
• minbatch: the @BATCHES batch size over the columns

Examples

julia> using SimilaritySearch

julia> X = rand(Float32, 128, 1000);

julia> hp = Projections.HadamardProjection(128);

julia> Y = Projections.transform(hp, X);

julia> size(Y)
(128, 1000)

transform(p::PCAProjection, v::AbstractVector)

Projects the vector v (of length indim(p)) onto p's principal directions, returning a new vector of length outdim(p).

transform(p::PCAProjection, X::AbstractMatrix)

Projects every column (vector) of X onto p's principal directions, returning a new matrix with outdim(p) rows and the same number of columns as X. Unlike RandomProjections/HadamardProjection, this is a single call into MultivariateStats (already a vectorized BLAS matrix-matrix product), so there is no minbatch to parallelize over.

Examples

julia> using SimilaritySearch

julia> X = rand(Float32, 128, 1000);

julia> p = SimilaritySearch.Projections.PCAProjection(X, 32);

julia> Y = SimilaritySearch.Projections.transform(p, X);

julia> size(Y)
(32, 1000)

transform!(rp::RandomProjections, out::AbstractVector, v::AbstractVector)

In-place version of transform: projects v using rp and stores the result in out, which must have length outdim(rp). Returns out.

Arguments

• rp: the projection to apply
• out: the output vector where the projected vector is stored
• v: the input vector to project, of length indim(rp)

transform!(rp::RandomProjections, O::AbstractMatrix, X::AbstractMatrix; minbatch::Int=4)

In-place version of transform(rp, X): projects every column of X using rp and stores the result in O, which must have outdim(rp) rows and the same number of columns as X. Returns O.

Arguments

• rp: the projection to apply
• O: the output matrix where the projected vectors are stored
• X: a matrix whose columns are the vectors to project, each of length indim(rp)
• minbatch: minimum number of columns processed per parallel task (see @BATCHES)

transform!(hp::HadamardProjection, out::AbstractVector, v::AbstractVector)

In-place version of transform: projects v using hp and stores the result in out, which must have length indim(hp) (== outdim(hp)). Returns out.

Arguments

• hp: the projection to apply
• out: the output vector where the projected vector is stored, of length indim(hp)
• v: the input vector to project, of length indim(hp)

transform!(hp::HadamardProjection, O::AbstractMatrix, X::AbstractMatrix; minbatch::Int=4)

In-place version of transform(hp, X): projects every column of X using hp and stores the result in O, which must have the same size as X. Returns O.

One fwht! per column, over @BATCHES of minbatch columns, the copy from X included: nothing is planned and nothing is locked, so the columns transform in parallel (a serial copyto!(O, X) ahead of the loop cost three times the transform itself). (Issue #54 had made this path a single batched FFTW call to avoid rebuilding a plan per column; issue #89 removed the plans altogether, and the per-column butterfly is 50x faster than that batched call.)

Arguments

• hp: the projection to apply
• O: the output matrix where the projected vectors are stored
• X: a matrix whose columns are the vectors to project, each of length indim(hp)
• minbatch: the @BATCHES batch size over the columns

transform!(p::PCAProjection, out::AbstractVector, v::AbstractVector)
transform!(p::PCAProjection, O::AbstractMatrix, X::AbstractMatrix)

In-place versions of transform: projects v/X using p and stores the result in out/O, which must have length/row-count outdim(p). Returns out/O.

bitsketch(R::AbstractMatrix{<:AbstractFloat}, v::AbstractVector{<:AbstractFloat}) -> Vector{UInt64}
bitsketch(R::AbstractMatrix{<:AbstractFloat}, X::AbstractMatrix{<:AbstractFloat}; minbatch::Int=4) -> Matrix{UInt64}

Computes a random-rotation bit sketch (a SimHash-style binary locality-sensitive hash): rotates the input by R (the same convention as RandomProjections, i.e. v -> R' * v) and encodes the sign of each of the size(R, 2) resulting coordinates as one bit – non-negative maps to 1, negative maps to 0 (see packsigns) – packed into UInt64 words (64 bits per word, cld(size(R, 2), 64) words per sketch). Vectors whose rotated coordinates fall on the same side of the random hyperplanes defined by the columns of R hash to the same bit pattern, so the Hamming distance between two sketches approximates the angular distance between the corresponding original vectors.

Arguments

• R: the rotation matrix, of size (indim, outdim); build one with gaussian or qr (pass its .map), or use bitsketch(method, outdim, data) to build and apply a fresh one in a single step
• v/X: the vector, or matrix (one vector per column), to sketch; must be given as Float32/Float64 (or another AbstractFloat subtype)
• minbatch: (matrix method only) minimum number of columns processed per parallel task

Examples

julia> using SimilaritySearch

julia> R = SimilaritySearch.Projections.gaussian(128, 256).map;

julia> X = rand(Float32, 128, 1000);

julia> B = SimilaritySearch.Projections.bitsketch(R, X);

julia> size(B), eltype(B)  # (4, 1000), UInt64  -- cld(256, 64) == 4 words per sketch

bitsketch(rp::RandomProjections, v::AbstractVector{<:AbstractFloat}) -> Vector{UInt64}
bitsketch(rp::RandomProjections, X::AbstractMatrix{<:AbstractFloat}; minbatch::Int=4) -> Matrix{UInt64}

Computes a bit sketch (see bitsketch(R::AbstractMatrix, data)) using an already-built RandomProjections rotation rp, equivalent to bitsketch(getmap(rp), data). Reuse the same rp to sketch a query the same way as an already-sketched dataset, so that the resulting sketches are comparable.

Examples

julia> using SimilaritySearch

julia> rp = SimilaritySearch.Projections.gaussian(128, 256);

julia> X = rand(Float32, 128, 1000);

julia> B = SimilaritySearch.Projections.bitsketch(rp, X);

julia> q = rand(Float32, 128);

julia> bq = SimilaritySearch.Projections.bitsketch(rp, q);  # same rotation, comparable to B's columns

bitsketch(method::Symbol, outdim::Int, data::AbstractVecOrMat{<:AbstractFloat};
          rng::AbstractRNG=Random.default_rng(), FloatType::Type=Float32, minbatch::Int=4)
    -> (bitsketch, R)

Convenience bitsketch that builds a fresh rotation matrix R – using gaussian (method = :gaussian) or qr (method = :qr) – applies it in a single step, and returns both the resulting sketch(es) and R as a tuple (bitsketch, R), so that the very same rotation can be reused afterwards (e.g. via bitsketch(R, data)) to sketch further vectors comparable to this one (e.g. queries against an already-sketched dataset).

Arguments

• method: :gaussian or :qr, selecting which of gaussian/qr builds the rotation matrix; any other value raises ArgumentError
• outdim: the number of sketch bits (i.e., the output dimension of the rotation)
• data: the vector or matrix (one vector per column) to sketch
• rng, FloatType: forwarded to the chosen rotation-matrix generator
• minbatch: (matrix method only) minimum number of columns processed per parallel task

Examples

julia> using SimilaritySearch

julia> X = rand(Float32, 128, 1000);

julia> B, R = SimilaritySearch.Projections.bitsketch(:gaussian, 256, X);

julia> size(B), eltype(B)  # (4, 1000), UInt64

julia> q = rand(Float32, 128);

julia> bq = SimilaritySearch.Projections.bitsketch(R, q);  # reuse R, comparable to B's columns

bitsketch(hp::HadamardProjection, v::AbstractVector{<:AbstractFloat}) -> Vector{UInt64}
bitsketch(hp::HadamardProjection, X::AbstractMatrix{<:AbstractFloat}; minbatch::Int=4) -> Matrix{UInt64}

Computes a bit sketch the same way as bitsketch(R::AbstractMatrix, data), but using the fast Walsh-Hadamard transform (HadamardProjection) instead of a dense random rotation matrix: it encodes the sign of each of the indim(hp) transformed coordinates into UInt64-packed bits (see packsigns).

Arguments

• hp: the HadamardProjection to apply
• v/X: the vector, or matrix (one vector per column), to sketch
• minbatch: (matrix method only) minimum number of columns processed per parallel task

Examples

julia> using SimilaritySearch

julia> hp = SimilaritySearch.Projections.HadamardProjection(128);  # 128 is a power of two

julia> X = rand(Float32, 128, 1000);

julia> B = SimilaritySearch.Projections.bitsketch(hp, X);

julia> size(B), eltype(B)  # (2, 1000), UInt64  -- cld(128, 64) == 2 words per sketch

bitsketch(p::PCAProjection, v::AbstractVector{<:AbstractFloat}) -> Vector{UInt64}
bitsketch(p::PCAProjection, X::AbstractMatrix{<:AbstractFloat}; minbatch::Int=4) -> Matrix{UInt64}

Computes a bit sketch the same way as bitsketch(R::AbstractMatrix, data), but using a fitted PCAProjection instead of a random rotation: it encodes the sign of each of the outdim(p) transformed coordinates into UInt64-packed bits (see packsigns). As with the random-rotation forms, reuse the same p to sketch a query and the dataset it will be compared against, so their sketches are comparable.

Arguments

• p: the PCAProjection to apply
• v/X: the vector, or matrix (one vector per column), to sketch
• minbatch: (matrix method only) minimum number of columns processed while packing sign bits (see packsigns); transforming X through p itself is not batched – see transform(p::PCAProjection, X::AbstractMatrix)

Examples

julia> using SimilaritySearch

julia> X = rand(Float32, 128, 1000);

julia> p = SimilaritySearch.Projections.PCAProjection(X, 256);

julia> B = SimilaritySearch.Projections.bitsketch(p, X);

julia> size(B), eltype(B)  # (4, 1000), UInt64  -- cld(256, 64) == 4 words per sketch

bitsketch(m::DistantHyperplanes, obj) -> Vector{UInt64}
bitsketch(m::DistantHyperplanes, X::AbstractDatabase; minbatch::Int=4) -> MatrixDatabase

Encodes obj (or every object of X) with the hyperplanes of m: bit i is 1 when obj falls on the "<=" side of hyperplane m.H[i] (see DistantHyperplanes), packed into UInt64 words. Sketches are compared with distance(m).

Arguments

• m: the fitted DistantHyperplanes sketch generator
• obj/X: the object, or database of objects, to sketch
• minbatch: (database method only) minimum number of items processed per parallel task

bitsketch(m::AnchoredDistantHyperplanes, obj) -> Vector{UInt64}
bitsketch(m::AnchoredDistantHyperplanes, X::AbstractDatabase; minbatch::Int=4) -> MatrixDatabase

Encodes obj (or every object of X) with the hyperplanes of m, exactly as bitsketch does for a DistantHyperplanes. Sketches are compared with distance(m).

Arguments

• m: the fitted AnchoredDistantHyperplanes sketch generator
• obj/X: the object, or database of objects, to sketch
• minbatch: (database method only) minimum number of items processed per parallel task

bitsketch(B::RandomHyperplanes, v) -> Vector{UInt64}
bitsketch(B::RandomHyperplanes, X::AbstractDatabase; minbatch::Int=4) -> MatrixDatabase

Encodes v (or every object of X) with the reference pairs of B (see RandomHyperplanes), packed into UInt64 words. Sketches are compared with distance(B).

Arguments

• B: the RandomHyperplanes sketch generator
• v/X: the object, or database of objects, to sketch
• minbatch: (database method only) minimum number of items processed per parallel task

packsigns(y::AbstractVector{<:Real}) -> Vector{UInt64}

Out-of-place version of packsigns!: packs the sign of each entry of y into a freshly allocated Vector{UInt64} of length cld(length(y), 64).

packsigns(Y::AbstractMatrix{<:Real}; minbatch::Int=4) -> Matrix{UInt64}

Packs the sign of every column of Y (see packsigns) into a Matrix{UInt64} with cld(size(Y, 1), 64) rows and the same number of columns as Y. Columns are packed in parallel using @BATCHES.

packsigns!(out::AbstractVector{UInt64}, y::AbstractVector{<:Real})

Packs the sign of each entry of y into out (a pre-allocated, cld(length(y), 64)-long UInt64 vector): entries >= 0 are encoded as bit 1, negative entries as bit 0. If length(y) is not a multiple of 64, the unused high bits of the last word are left as 0. Returns out.

HadamardProjection(indim::Int)
HadamardProjection(indim::Int, outdim::Int)

Wraps a fast Walsh-Hadamard transform (FWHT), used as an orthogonal change of basis (via transform/transform!), analogous in purpose to RandomProjections but computed with the $O(n \log n)$ FWHT (an in-place butterfly, fwht!; no plan, no FFTW) instead of a dense matrix-vector product, and requiring no random matrix to be generated or stored.

Unlike RandomProjections, HadamardProjection does not reduce dimensionality: transform always returns as many coordinates as it received (outdim(hp) == indim(hp)), since fwht_natural! computes a full, exact (up to normalization), orthogonal transform of its input, in natural Hadamard ordering. The two-argument constructor exists only to make outdim explicit at call sites that already pass one to other projection types (e.g. RandomProjections); it requires outdim == indim and raises ArgumentError otherwise.

Arguments

• indim: the dimension of the input vectors; must be a power of two (fwht_natural! requirement), otherwise an ArgumentError is thrown
• outdim: if given, must equal indim (otherwise an ArgumentError is thrown), since this projection does not support dimensionality reduction/truncation

Examples

julia> using SimilaritySearch

julia> hp = Projections.HadamardProjection(128);

julia> Projections.indim(hp), Projections.outdim(hp)
(128, 128)

fwht!(a::AbstractVector) -> a

The fast Walsh-Hadamard transform of a in place, in natural (Hadamard) ordering and scaled by 1/length(a): the same values Hadamard.jl's fwht_natural! produces, bit for bit, which is what every encoding built on HadamardProjection was measured with. length(a) must be a power of two (unchecked here; the projection's constructor checks).

It is the plain iterative butterfly, log2(n) passes of sums and differences over halves of growing size, and it needs no plan. That is the point of it: the transform used to go through FFTW, which built a plan under its global planning lock on every per-vector call (170-400 µs per vector at 128-4096 dimensions, and worse with threads, which serialize on the lock), and whose execution of this shape – a 2×2×…×2 multidimensional r2r transform – stays at 3-54 µs even with a cached, MEASURE-quality plan (issue #89). Don't bring a planned transform back for this.

For a contiguous Float32/Float64 vector (a Vector, or a column view of a matrix) the first three passes, whose inner loops have 1, 2 and 4 elements and cannot vectorize, are one sweep of 8-point blocks done with three lane shuffles of a Vec{8}, the remaining passes load and store Vec{8}s, and the 1/n scale rides in the last pass: 95 ns at 128 dimensions, 0.4 µs at 512, 1.0 µs at 1024 and 4.6 µs at 4096 in Float32 (33-800 ns per column of a matrix over 64 threads), 4-5x the scalar loops, which any other element type or layout still gets. Every variant is bit-identical: each output is the same sum or difference of the same two values, in the same order.

RandomizedHadamard(dim; rng=Random.default_rng())

A random orthogonal rotation computed in dim log dim flops: a random sign flip of every coordinate followed by the Walsh-Hadamard transform (Projections.HadamardProjection, which scales by 1/dim), rescaled by sqrt(dim) so that norms are preserved. dim must be a power of two. The plain transform is not a random rotation (a constant vector lands on a single coordinate), and the sign flip is what makes it one.

Rotation

What the estimators take as their rotation: a RandomProjections with a square map, such as Projections.qr(dim, dim) (an orthogonal matrix, dim² floats and dim² flops per vector), or a RandomizedHadamard (dim a power of two, dim signs stored, dim log dim flops). Both implement rotate(rot, v) and rotationdim(rot); ScalarQuant.SQEncoder also takes nothing, which rotates nothing.

PCAProjection(X::AbstractMatrix{<:AbstractFloat}, outdim::Int; kwargs...)

Wraps a PCA-based dimensionality reduction (via MultivariateStats.jl's PCA), analogous in purpose to RandomProjections but fitted from data instead of drawn at random: transform/transform! projects a vector/matrix onto the outdim orthogonal directions of largest variance in X (after centering by X's mean). Unlike a random projection, the reduction this produces depends on, and is tailored to, the specific data it was fitted on – so, unlike RandomProjections/ HadamardProjection, it cannot be regenerated independently of that data; reuse the same PCAProjection (as with those, e.g. via transform) to project a query comparably to an already-projected dataset.

Arguments

• X: the training data, one object per column, used to fit the principal directions
• outdim: the number of output dimensions to keep (MultivariateStats' maxoutdim)
• kwargs...: forwarded to MultivariateStats.fit(MultivariateStats.PCA, X; maxoutdim=outdim, kwargs...) (e.g. method=:svd/:cov, pratio, mean)

Examples

julia> using SimilaritySearch

julia> X = rand(Float32, 128, 1000);

julia> p = SimilaritySearch.Projections.PCAProjection(X, 32);

julia> indim(p), outdim(p)
(128, 32)

julia> Y = SimilaritySearch.Projections.transform(p, X);

julia> size(Y)
(32, 1000)

DistantHyperplanes(dist::SemiMetric, X::AbstractDatabase, nbits::Int;
    hsel::Int=nbits*1024, entquantile::Float64=0.9, henc::Int=2^13,
    minbatch::Int=4, verbose::Bool=true)

Builds a hyperplane-based binary sketch generator for a generic metric space. A hyperplane is a pair of anchor objects (i, j) sampled from X; an object obj falls on one side of hyperplane (i, j) when evaluate(dist, obj, X[i]) <= evaluate(dist, obj, X[j]).

hsel candidate hyperplanes are sampled and characterized by which side a random subsample of henc objects from X falls on; hyperplanes whose resulting bit-vector is not close to a fair coin (entropy, out of a maximum of 1.0) are discarded as uninformative. The entropy cutoff is the entquantile quantile of the entropies the candidates actually achieved, i.e. the top 1 - entquantile fraction of candidates (by entropy) survive – this adapts automatically to whatever entropy ceiling the dataset/distance can actually deliver, instead of comparing against a fixed absolute bar that could sit above that ceiling and silently discard every candidate, producing a useless 0-bit sketch. Only a fully degenerate case – every candidate has exactly zero entropy – raises an error, since no cutoff can rescue it. From the survivors, nbits hyperplanes are finally kept, chosen to be as mutually diverse as possible – i.e., spread apart in Hamming distance between their characterization bit-vectors, up to a global bit-flip (two hyperplanes whose sides happen to be labeled oppositely are just as good a pair to keep as two that agree) – via a farthest-first traversal (fft).

Use bitsketch to encode objects with the resulting DistantHyperplanes; the sketches it produces are compared with distance(m) (Hamming distance over UInt64-packed bits).

Arguments

• dist: the distance function of the underlying metric space
• X: the database used both to sample candidate hyperplane anchors and, later, as the reference set hyperplanes are evaluated against when encoding new objects
• nbits: the number of output bits (i.e., hyperplanes) to keep; must be a multiple of 64

Keyword Arguments

• hsel: number of candidate hyperplanes (pairs of objects) to sample and characterize
• entquantile: quantile (in [0, 1]) of the candidates' achieved entropies used as the acceptance cutoff (see above)
• henc: sample size used to characterize each candidate hyperplane; must be a multiple of 64 and smaller than length(X)
• minbatch: minimum number of items processed per parallel task (see @BATCHES)
• verbose: whether the per-center progress message of the underlying fft call is produced

Examples

julia> using SimilaritySearch

julia> X = MatrixDatabase(rand(Float32, 8, 10_000));

julia> m = SimilaritySearch.Projections.DistantHyperplanes(SimilaritySearch.Dist.L2(), X, 128);

julia> B = SimilaritySearch.Projections.bitsketch(m, X);

julia> size(B.matrix), eltype(B.matrix)  # (2, 10000), UInt64 -- 128 bits / 64 = 2 words per sketch

AnchoredDistantHyperplanes(dist::SemiMetric, X::AbstractDatabase, nbits::Int;
    anchor=nothing, anchorpolicy::Symbol=:random,
    hsel::Int=nbits*1024, entquantile::Float64=0.9, henc::Int=2^13,
    minbatch::Int=4, verbose::Bool=true)

A variant of DistantHyperplanes that avoids the flip ambiguity of hyperplanes by anchoring, instead of masking it after the fact.

A hyperplane defined by anchors (i, j) is, geometrically, the very same hyperplane as (j, i): relabeling which point is "first" just flips every bit of its characterization (the side called 1 becomes the side called 0). Plain candidate sampling has no reason to prefer one labeling over the other, so DistantHyperplanes has to compare candidates with a flip-invariant distance (min(hamming(u, v), hamming(u, ~v))) to avoid mistaking two near-identical hyperplanes – that only happen to disagree on which side got called 1 – for a genuinely diverse pair.

AnchoredDistantHyperplanes sidesteps this by fixing the labeling convention up front: a single reference object, the anchor, orders every candidate pair (i, j) so that i is always the one closer to the anchor (swapping them if sampling produced the opposite order). With that shared convention in place, "closer to i" consistently means "closer to the anchor-proximal point" for every hyperplane, so plain Hamming distance – not the flip-invariant one – is enough to tell diverse hyperplanes from redundant ones.

Whether a given anchor choice helps or hurts is an open, dataset-dependent question – anchor accepts an explicit choice (an object, or an integer id into X) for that reason; when left unset, one is picked automatically per anchorpolicy:

• :random (default): a uniformly random object of X
• :extremal: the farthest object, in dist, from a random starting point (one step of fft) – tends to sit at the periphery of X, which may spread distances-to-anchor out more than a typical (e.g. random) point would

Use bitsketch to encode objects, exactly as with DistantHyperplanes; sketches are compared with distance(m).

Arguments

• dist: the distance function of the underlying metric space
• X: the database used both to sample candidate hyperplane anchors and, later, as the reference set hyperplanes are evaluated against when encoding new objects
• nbits: the number of output bits (i.e., hyperplanes) to keep; must be a multiple of 64

Keyword Arguments

• anchor: the anchor object to orient hyperplane pairs by; an integer is taken as an id into X (i.e., X[anchor]), anything else is used directly as the anchor object. nothing (default) computes one automatically, per anchorpolicy
• anchorpolicy: :random or :extremal (see above); only used when anchor === nothing
• hsel: number of candidate hyperplanes (pairs of objects) to sample and characterize
• entquantile: quantile (in [0, 1]) of the candidates' achieved entropies used as the acceptance cutoff – adapts automatically to whatever entropy ceiling the dataset/distance can actually deliver, instead of comparing against a fixed absolute bar that could silently discard every candidate (see DistantHyperplanes and issue #55)
• henc: sample size used to characterize each candidate hyperplane; must be a multiple of 64 and smaller than length(X)
• minbatch: minimum number of items processed per parallel task (see @BATCHES)
• verbose: whether the per-center progress message of the underlying fft call is produced

Examples

julia> using SimilaritySearch

julia> X = MatrixDatabase(rand(Float32, 8, 10_000));

julia> m = SimilaritySearch.Projections.AnchoredDistantHyperplanes(SimilaritySearch.Dist.L2(), X, 128);

julia> m2 = SimilaritySearch.Projections.AnchoredDistantHyperplanes(SimilaritySearch.Dist.L2(), X, 128; anchor=1); # X[1] as anchor

julia> B = SimilaritySearch.Projections.bitsketch(m, X);

RandomHyperplanes(dist::SemiMetric, refs::AbstractDatabase, npairs::Integer)

Binary sketch generator for a generic metric space, based on npairs pairs of reference objects (anchors) drawn from refs – so refs must hold exactly 2npairs objects, laid out as consecutive pairs (refs[2i-1], refs[2i] is the i-th pair). An object is encoded into npairs bits: bit i is 1 when the object is closer to refs[2i-1] than to refs[2i] (see bitsketch).

Unlike DistantHyperplanes, which searches for and filters informative, mutually-diverse anchor pairs from data, RandomHyperplanes takes the anchor pairs as given – e.g., a plain random sample of the dataset – making it a much cheaper, simpler baseline.

Sketches are compared with distance(B) (Hamming distance over UInt64-packed bits).

Arguments

• dist: the distance function of the underlying metric space
• refs: the 2npairs reference objects, laid out as consecutive pairs
• npairs: the number of output bits (i.e., reference pairs); must be a multiple of 64

Examples

julia> using SimilaritySearch

julia> X = MatrixDatabase(rand(Float32, 8, 10_000));

julia> refs = SubDatabase(X, rand(1:10_000, 256));  # 128 pairs -> 128 bits

julia> m = SimilaritySearch.Projections.RandomHyperplanes(SimilaritySearch.Dist.L2(), refs, 128);

julia> B = SimilaritySearch.Projections.bitsketch(m, X);

julia> size(B.matrix), eltype(B.matrix)  # (2, 10000), UInt64 -- 128 bits / 64 = 2 words per sketch

QuantSketch{NBITS,MODEL}

An NBITS-per-component sketch encoder: it wraps any sketch model that sketchvalues! understands and, instead of keeping one sign bit per component the way bitsketch does, keeps an NBITS-wide unsigned code per component, quantized with the global (database-wide) scalar quantizers SQgu2/SQgu4/SQgu8.

The point is that a sign bit throws away how far an object sits from each hyperplane: an object hugging a hyperplane and one far across it get the same bit, so Hamming distance cannot tell a marginal disagreement from a decisive one. An NBITS code keeps that magnitude, and the sketches are then compared with the squared-Euclidean distance over the packed codes (distance(qs)), which the SQgu* kernels evaluate with SIMD directly on the codes – no dequantization – because a single global affine map is shared by every component of every sketch.

NBITS == 1 is supported and is exactly bitsketch's encoding (UInt64-packed sign bits, compared with Hamming), so a sweep over 1, 2, 4, 8 bits runs through one uniform API instead of two.

Encode with quantsketch(qs, obj_or_db), compare with distance(qs), and see SketchedSearch for the ready-made index-and-query pipeline built on top.

Memory

At m = outdim(model) components a sketch occupies m * NBITS / 8 bytes, i.e. NBITS times what the binary sketch of the same model costs. A fair comparison across widths therefore holds bits fixed, not components: 1024 binary hyperplanes, 512 at 2 bits, 256 at 4 bits and 128 at 8 bits all occupy 128 bytes.

QuantSketch(model, nbits::Int, X;
            minmax=nothing, normalize::Bool=true,
            quant=[0.025, 0.975], samplesize::Int=0)

Fits an encoder for model at nbits bits per component, estimating the global quantization range from X.

Arguments

• model: the sketch model; anything sketchvalues! accepts
• nbits: 1, 2, 4 or 8; anything else raises ArgumentError
• X: the dataset the range is estimated from – an AbstractDatabase or a matrix whose columns are the objects. Ignored when nbits == 1 (a sign needs no range) or when minmax is given explicitly.

Keyword Arguments

• minmax: an explicit (min, max) global range, skipping the estimation from X
• normalize: whether to divide each component by its hyperplane's width (see hyperplanewidths) before quantizing, putting every component on a common [-1, 1] scale. Defaults to true and is a no-op for models with no such width (the rotation/projection family).
• quant: the lower/upper quantiles of the sampled values used to estimate minmax; quantiles rather than extrema so a few outliers cannot flatten the useful range
• samplesize: how many objects of X are characterized to estimate minmax; 0 (the default) picks clamp(ceil(Int, sqrt(length(X))), 1, length(X))

Examples

julia> using SimilaritySearch, SimilaritySearch.Projections

julia> X = MatrixDatabase(rand(Float32, 8, 10_000));

julia> m = DistantHyperplanes(SimilaritySearch.Dist.L2(), X, 128; verbose=false);

julia> qs = QuantSketch(m, 4, X);          # 4 bits per hyperplane instead of 1

julia> B = quantsketch(qs, X);             # MatrixDatabase of packed codes

julia> distance(qs)                        # SQgu4.SqL2(), evaluated on the packed codes

quantsketch(qs::QuantSketch, obj) -> Vector
quantsketch(qs::QuantSketch, X::AbstractDatabase; minbatch::Int=4) -> MatrixDatabase
quantsketch(qs::QuantSketch, X::AbstractMatrix; minbatch::Int=4) -> Matrix

Encodes obj (or every object of X) with qs, returning sketchsize(qs)-long codes – UInt64 sign words at 1 bit, packed UInt8 codes otherwise. Sketches are comparable only against others produced by the same qs (the quantization range and, when normalize is on, the per-hyperplane scales are baked into it), and must be compared with distance(qs).

The collection methods mirror bitsketch's return conventions – a database in, a MatrixDatabase out; a matrix in, a Matrix out – and never materialize the dense (outdim, n) matrix of real values: each batch keeps a single outdim-long scratch vector reused across the objects it encodes.

Arguments

• qs: the fitted encoder
• obj/X: the object, or collection of objects, to sketch
• minbatch: (collection methods only) minimum number of objects processed per parallel task

quantsketch(model, nbits::Int, X; minbatch::Int=4, kwargs...) -> (codes, qs)

Convenience one-step form: fits a QuantSketch for model at nbits bits over X and immediately encodes X with it, returning both the codes and the encoder, so the very same encoder can sketch further objects (queries, above all) comparably. Extra keyword arguments are forwarded to QuantSketch.

Examples

julia> using SimilaritySearch, SimilaritySearch.Projections

julia> X = MatrixDatabase(rand(Float32, 8, 10_000));

julia> m = RandomHyperplanes(SimilaritySearch.Dist.L2(), SubDatabase(X, rand(1:10_000, 256)), 128);

julia> B, qs = quantsketch(m, 4, X);

julia> bq = quantsketch(qs, X[1]);   # a query, encoded the same way

sketchvalues!(out::AbstractVector{Float32}, model, obj) -> out
sketchvalues(model, obj) -> Vector{Float32}

The real-valued vector a sketch model produces for obj before it is discretized – the quantity bitsketch reduces to one sign bit per component and QuantSketch keeps at 2, 4 or 8 bits instead. out must have length outdim(model).

Two families of models are supported, with a deliberately shared sign convention (a component is >= 0 exactly when bitsketch would set its bit to 1):

• rotations/projections (RandomProjections, HadamardProjection, PCAProjection): the projected coordinates themselves, i.e. transform!. The magnitude says how far the object sits from the hyperplane normal to that direction.
• metric hyperplanes (DistantHyperplanes, AnchoredDistantHyperplanes, RandomHyperplanes): the signed margin of obj against each hyperplane (a, b), d(obj, b) - d(obj, a), which is positive when obj is closer to a (the side those models encode as bit 1) and whose magnitude says how decisively so.

Arguments

• out: the output vector, of length outdim(model)
• model: the projection or hyperplane model
• obj: the object to characterize

sketchbits(qs::QuantSketch)

The number of bits qs spends on each component (1, 2, 4 or 8).

sketchsize(qs::QuantSketch)

The length, in storage words, of one sketch produced by qs: cld(outdim(qs), 64) UInt64 words at 1 bit, and cld(outdim(qs), 8 ÷ nbits) UInt8 bytes otherwise.

hyperplanewidths(model) -> Vector{Float32}

The separation d(a, b) between the two anchors of every hyperplane of a metric hyperplane model (DistantHyperplanes, AnchoredDistantHyperplanes, RandomHyperplanes). By the triangle inequality the margin of any object against hyperplane i (see sketchvalues!) lies in [-w[i], w[i]], so dividing by w puts every component of the sketch on the same, dimensionless [-1, 1] scale – which is what makes a single global quantization range (the SQgu2/SQgu4/SQgu8 family, and hence their SIMD code-space distances) meaningful across hyperplanes built from anchor pairs that are, individually, arbitrarily far apart or close together. See QuantSketch's normalize keyword.

Returns nothing for models that have no such natural per-component width (the rotation/projection family, whose coordinates already share one scale).

SketchedSearch{QS,IDX,DIST,DB} <: AbstractSearchIndex

SketchedSearch(model, nbits::Int, dist::PreMetric, db::AbstractDatabase;
               factor::Int=8, index=exhaustivesketchindex, minbatch::Int=4, kwargs...)

The whole sketch-based search pipeline behind one ordinary index: it sketches db with model at nbits bits per component (see QuantSketch), indexes the sketches, and answers a query in the original space by sketching it the same way, collecting factor * k candidates under the cheap sketch distance, and re-scoring exactly those with dist.

Assembling that by hand is only a handful of lines, but every one of them is a place to get it wrong silently: sketching the query with a freshly-fitted encoder instead of the dataset's (incomparable codes), indexing the sketches under the original distance instead of distance(qs), or forgetting the re-scoring pass and returning the sketch distances as if they were real ones. None of those raise an error – they just quietly return worse results. Going through SketchedSearch makes them unrepresentable.

SketchedSearch is built once over a fixed db and is not incremental: its encoder's quantization range (and, for the hyperplane models, its anchors) are fitted on exactly the data given here, so growing the collection means rebuilding it rather than pushing into it.

Because it is an AbstractSearchIndex, search/searchbatch work on it unchanged, and the ids and distances it returns are ids into db and true dist values, so it is a drop-in replacement for an exact index – swapping one in is the whole experiment. factor is the knob that trades recall for time: the sketch stage is fast but lossy, so it is asked for more candidates than needed and the exact distance settles the final order among them.

Arguments

• model: the sketch model; anything sketchvalues! accepts (a hyperplane model like DistantHyperplanes, or a rotation like RandomProjections)
• nbits: bits per component – 1 (plain bitsketch + Hamming), 2, 4 or 8
• dist: the original, exact distance, used to re-score candidates
• db: the original database

Keyword Arguments

• factor: how many times k candidates the sketch stage retrieves before re-scoring; 1 disables the widening (but still re-scores, so distances stay exact)
• index: (sketchdist, sketches) -> AbstractSearchIndex, building the index over the sketches. Defaults to an ExhaustiveSearch – a linear scan, which is the point: it is over sketches many times smaller than the objects. Pass a closure returning a ready-to-search index (already index!-ed) to use e.g. a SearchGraph instead.
• minbatch: minimum number of objects per parallel task while sketching db
• remaining keywords are forwarded to QuantSketch (normalize, minmax, quant, samplesize)

Examples

julia> using SimilaritySearch, SimilaritySearch.Projections

julia> dist = SimilaritySearch.Dist.L2();

julia> db = MatrixDatabase(rand(Float32, 8, 10_000));

julia> Q  = MatrixDatabase(rand(Float32, 8, 100));

julia> m = DistantHyperplanes(dist, db, 128; verbose=false);

julia> S = SketchedSearch(m, 4, dist, db);          # 4-bit codes, 8x candidate widening

julia> ids, dists = searchbatch(S, GenericContext(), Q, 10);   # ids into db, exact dists

exhaustivesketchindex(sketchdist, sketches) -> ExhaustiveSearch

The default index builder of SketchedSearch: a plain ExhaustiveSearch over the sketches. A linear scan is a sensible default here precisely because the sketches are far smaller than the objects they stand for, and it keeps the recall of the pipeline attributable to the sketch alone – no second layer of approximation on top.