Verifying NTQR¶
NTQR verifies evaluations, who or what verified NTQR? Althought the mathematics of computing the consistent evaluation set has only linear algebra operations that are relatively easy to understand, the implementation of that math may be buggy or incorrect.
There are four properties that we would any implementation to obey:
Validity
Completeness
Uniformity
Mixing
Hampering validation of 2,3, and 4 is the large sizes of the evaluation sets involved. To understand how large these sets can get let us look at the corner case of only one label in the answer key. This is single label state for the answer key is represented in the geometry of NTQR as vertices on a simplex.
Validity: Joint evaluations must obey axioms¶
Validity means that points returned by any of NTQR generators must be valid, that they obey the axioms of their set. For the PossibleSet that means that label responses must sum to the count of the label in the answer key, the label simplex axioms. Evaluations in the ConsistentSet must obey the \(R\) label axioms and the test event axioms (sum of event counts across labels equals event test count). Validity is thus a set of linear sums and easy to perform.
Let us do so.
import itertools, random, collections
import sympy
import ntqr
import ntqr.raxioms
import ntqr.evaluations
# We create synthetic test results for N=4 classifiers
# doing R=5 classification
labels = ('a', 'b', 'c','d','e')
classifiers = (1,2,3,4)
# The test will have Q=100 questions. NTQR
# is a logic so we randomly pick accuracy
# on our test sampling code to be uniform
Q=100
answer_key = [random.choice(labels) for i in range(Q)]
label_accuracies = [{label:random.uniform(0.0,1.0) for label in labels} for classifier in classifiers]
test_results = [tuple(tlabel if label_accuracies[j][tlabel]>random.random()
else random.choice([ol for ol in labels if ol != tlabel])
for j in range(len(classifiers)) )
for tlabel in answer_key]
counts=collections.Counter(test_results)
print(len(counts))
cSet=ntqr.evaluations.ConsistentSet(labels,classifiers,counts)
# And now we want to look at two sizes to make relative comparisons
# about how the position in the Q-simplex affects the evaluation set counts.
# We are looking at the PossibleSet which is huge. But luckily, because
# it is a uniform figure in response space, it is easy to compute its
# member size.
# We are using this fast computation of the possible set to 'probe'
# the density of evaluation sets in the Q-simplex
pSet = ntqr.evaluations.PossibleSet(labels,classifiers)
qCenter = (20,20,20,20,20)
qbVertex = (0,100,0,0,0)
print(f"Possible set at Q-simplex centroid: {float(pSet.set_count_at_ql(qCenter))}")
print(f"Possible set at Q-simplex vertex: {float(pSet.set_count_at_ql(qbVertex))}")
cSet = ntqr.evaluations.ConsistentSet(labels,classifiers,counts)
Possible set at Q-simplex centroid: 2.0452481236353108e+188
Possible set at Q-simplex vertex: 7.73312619258231e+124
# The vertices in the Q-simplex also have the property
# that ConsistentSets are small around them, relatively
# speaking. In fact, there can only be one member of the
# consistent set at each vertex of the Q-simplex. Lets
# confirm it.
eval_points=list(eval_point for eval_point in cSet.set_generator(qbVertex))
eval_points
[(HashablePoint(<Compressed Sparse Row sparse matrix of dtype 'int64'
with 0 stored elements and shape (1, 625)>),
HashablePoint(<Compressed Sparse Row sparse matrix of dtype 'int64'
with 67 stored elements and shape (1, 625)>),
HashablePoint(<Compressed Sparse Row sparse matrix of dtype 'int64'
with 0 stored elements and shape (1, 625)>),
HashablePoint(<Compressed Sparse Row sparse matrix of dtype 'int64'
with 0 stored elements and shape (1, 625)>),
HashablePoint(<Compressed Sparse Row sparse matrix of dtype 'int64'
with 0 stored elements and shape (1, 625)>))]
# Is this point a valid member?
# Let's check PossibleSet membership
pSet.is_valid_point(eval_points[0],qbVertex)
(True, 'Valid')
# This point should also belong to the ConsistentSet
cSet.is_valid_point(eval_points[0], qbVertex)
(True, 'Valid')
Testing completeness of the generators¶
However large, the PossibleSet and the ConsistentSet are finite objects. Always. Do the generators in NTQR return all members? Are they complete?
The PossibleSet.set_generator code relies on code that produces all the points in a single label simplex. If that code is correct, barring implementation bug, the multi-label generator is complete. Furthermore, the high symmetry of PossibleSet means we have the function, PossibleSet.set_count_at_ql, that has a closed-form function that can compute the size of the set. So we can check for completeness of PossibleSet generators with these expected counts.
Checking completeness for ConsistentSet generators is harder. The observed counts break the high symmetry of PossibleSet and no closed form count of the set is possible. Here we must rely on corner cases that we can reason about easily.
These ‘corner cases’ are exactly that, around the vertices of the Q-simplex the ConsistentSet is small. At each vertex, there can only be one member in the set – the vector of the observed counts all assigned to a single label. One point off that is also easy to reason about. There should be \(E\) evaluations, where \(E\) are the observed events. Let’s check this.
# Let's count the sets by exhausting the generator
c_points = list(point for point in cSet.set_generator((1,99,0,0,0)))
f"Number of distinct decision events: {len(counts)}, set size: {len(c_points)}"
'Number of distinct decision events: 67, set size: 67'
Testing completeness of the random generators for ConsistentSet is also harder. Outside the simple cases near the Q-simplex vertex, we can only compare the deterministic and random generators.
Testing uniformity of the random set generators¶
import numpy as np
from scipy.stats import chisquare
import sys
def get_chain(generator, n_samples=100):
"""Get list."""
return list(itertools.islice(generator,n_samples))
def test_set_uniformity(samples):
"""
Checks if set members are returned uniformly.
Args:
samples: A list of evaluation points
"""
# 1. Count occurrences of each unique sample
counts = {}
for point in samples:
counts.setdefault(point,0)
counts[point] += 1
unique_points = len(counts)
# 2. Define observed and expected frequencies
# If the sampler is uniform, every possible state should be visited
# roughly the same number of times (total_samples / num_possible_states)
f_obs = [count for count in counts.values()]
total_samples = len(samples)
f_exp = np.full(unique_points, total_samples / unique_points)
print(unique_points)
# 3. Perform the Pearson's Chi-Squared test
statistic, p_value = chisquare(f_obs, f_exp=f_exp)
# Interpretation:
# p_value > 0.05: Fail to reject the null hypothesis;
# the distribution is consistent with uniformity.
# p_value < 0.05: Reject the null hypothesis;
# the samples show significant non-uniformity.
return statistic, p_value
points = get_chain(cSet.random_set_generator((1,1,97,1,0)),20_000)
len(points[200])
5
test_set_uniformity(points)
19343
(np.float64(646.7182), np.float64(1.0))
Testing mixing¶
import numpy as np
from sklearn.decomposition import PCA
def to_flat_array(obj):
# If it's a sparse matrix, convert to dense first
if hasattr(obj, "toarray"):
return obj.toarray().ravel()
# Otherwise, assume it's a standard numpy-like array
return np.array(obj).ravel()
def compute_r_hat_pca(chains_samples: list[list[tuple]]) -> float:
"""
Performs PCA-based convergence check on a product space.
Args:
chains_samples: A list of M chains, where each chain is a list
of N samples (each sample is a tuple of arrays).
Returns:
r_hat: The Gelman-Rubin statistic. Values close to 1.0
indicate convergence.
"""
# 1. Flatten all tuples in all chains to create a high-dimensional vector per sample
# Shape of all_vectors: (M * N, total_dimensions)
all_vectors = []
for chain in chains_samples:
for sample in chain:
# Concatenate all arrays in the tuple into one long flat vector
all_vectors.append(np.concatenate([to_flat_array(arr) for arr in sample]))
all_vectors = np.array(all_vectors)
# 2. Project into 1D space using PCA
# This identifies the axis of maximum variance across the entire pooled dataset
pca = PCA(n_components=1)
projected = pca.fit_transform(all_vectors)
# 3. Reshape back into (M chains, N samples per chain)
M = len(chains_samples)
N = len(chains_samples[0])
data = projected.reshape(M, N)
# 4. Calculate Gelman-Rubin (R-hat)
# B = Between-chain variance, W = Within-chain variance
chain_means = np.mean(data, axis=1)
grand_mean = np.mean(chain_means)
# B is the variance of chain means multiplied by N
B = (N / (M - 1)) * np.sum((chain_means - grand_mean)**2)
# W is the mean of individual chain variances
W = (1 / M) * np.sum(np.var(data, axis=1, ddof=1))
# var_theta is the estimated marginal posterior variance
var_theta = ((N - 1) / N) * W + (1 / N) * B
r_hat = np.sqrt(var_theta / W)
return r_hat
def get_chain(generator, n_samples=10):
"""Get list."""
return list(itertools.islice(generator,n_samples))
The mixing property is easier to check since it is capturing the geometry of the sets involved by using PCA. So, whereas in empirical checks of uniformity where we must deal with consistent sets small enough that we can realistically sample many times their size, we can try any Q-simplex point we want to check.
Let’s try the Q-simplex centroid.
# The mixing property is easier to check
qCentroid = (20,20,20,20,20)
chains = [get_chain(cSet.random_set_generator(qCentroid),10) for _ in range(100)]
compute_r_hat_pca(chains)
np.float64(1.0594061858211088)
chains = [get_chain(cSet.random_set_generator(qCentroid),100) for _ in range(100)]
compute_r_hat_pca(chains)
np.float64(1.0174249772909854)
This is pretty good since values of less than 1.1 are considered the standard for the Gelman-Rubin test in statistical software. However, these are new tests and code for me, so I may be using the test improperly. For that reason, none of the code in this notebook is in NTQR itself.
The development path for verifying NTQR¶
This notebook should make clear that the verification of NTQR is patchy and needs more work. My process is informal. The verification has started in this notebook in a manner that is easy to understand and hopefully will elicit user feedback on bugs, improvements, and alternatives.
It gives me a great thrill to see, in specific detail, what formal verification entails. The logic that NTQR is trying to implement correctly is simple to understand in its axioms and computations. How that logic is implemented to make it practical has, itself, raised verification issues. The who-judges-the-judges? problem is inescapable.