Price index methods are being used or currently being implemented by many statistical agencies around the world (e.g to compute the Consumer Price Index (CPI)) using bilateral or multilateral methods. For years the UK CPI was based on a bilateral method, with alternative data sources such as scanner or web-scraped data waiting in the wings. That transition has now happened: the ONS introduced groceries scanner data into consumer price inflation statistics in March 2026 (first incorporated from the February 2026 index), following the groceries scanner data analysis published in April 2025 and the overview of how scanner data are used published in January 2026. Rail fares and second-hand cars had already moved to alternative sources earlier in the transformation programme. The methodology the ONS settled on is GEKS-Törnqvist (internationally also known as CCDI) on a 25-month rolling window with a mean splice — not Geary-Khamis. Some context on why, and on where GK still matters, is at the end of this post.

Multilateral methods simultaneously make use of all data over a given time period. The use of multilateral methods for calculating temporal price indices is relatively new internationally, but these methods have been shown to have some desirable properties relative to their bilateral method counterparts, in that they account for new and disappearing products (to remain representative of the market) while also reducing the scale of chain-drift.

The Geary-Khamis method is a multilateral method which involves calculating a set of quality adjustment factors, $b_n$, simultaneously with the price levels, $P_t$. The two equations that determine both of these are:

\[\begin{aligned} b_{n}&=\sum_{t=1}^{T}\left[\frac{q_{t n}}{q_{n}}\right]\left[\frac{p_{t n}}{P_{t}}\right] \nonumber \\ P_{t}&=\frac{p^{t} \cdot q^{t}}{ \vec{b} \cdot q^{t}} \end{aligned}\]

These equations can be solved by an iterative method, where a set of $b_n$ are arbitrarily chosen, which can then be used to calculate an initial vector of price levels. This vector of prices is then used to generate a new $b$ vector, and so on until the changes become smaller than some threshold.

An alternative method is to use matrix operations, which is much more efficient. The problem of finding the $\vec{b}$ vector can be solved using the following system of equations:

\[\left[I_{N}-C\right] \vec{b}=0_{N}\]

where $I_n$ is the $N \times N$ identity matrix, $0_N$ is an $N$-dimensional vector of zeroes and $C$ is the $N \times N$ matrix given by,

\[C=\hat{q}^{-1} \sum_{t=1}^{T} s^{t} q^{t \mathbf{T}}\]

where $\hat{q}^{-1}$ is the inverse of an $N \times N$ diagonal matrix $\hat{q}$, where the diagonal elements are the total quantities purchased for each good over all time periods, $s^{t}$ is a vector of the expenditure shares for time period $t$, and $q^{t \mathbf{T}}$ is the transpose of the vector of quantities purchased in time period $t$. It can be shown that the matrix $[I - C]$ is singular since:

\[\bigg(\sum_{t=1}^T q^{t \mathbf{T}}\bigg) C = \sum_{t=1}^T q^{t \mathbf{T}}\]

So a normalisation is required to solve for $\vec{b}$ and we can use the constraint:

\[\sum_{n=1}^{N} b_{n} q_{n}=1\]

which is, in matrix form can be expressed as,

\[\vec{c}=R\left[\begin{array}{c} b_{1} q_{1} \\ \vdots \\ b_{N} q_{N} \end{array}\right] = \left[\begin{array}{c} 1 \\ 0 \\ \vdots \\ 0 \end{array}\right]\]

where $\vec{c}$ is an $N \times 1$ vector and the $R$ is the $N \times N$ matrix,

\[R=\left[\begin{array}{cccc} 1 & 1 & \ldots & 1 \\ 0 & \ldots & \ldots & 0 \\ \vdots & & & \vdots \\ 0 & \ldots & \ldots & 0 \end{array}\right]\]

Adding the constraint to the original equation we now have the solution for $\vec{b}$,

\[\vec{b}=\left[I_{N}-C+R\right]^{-1} \vec{c}\]

Once the $\vec{b}$ vector has been calculated, the price levels can be computed from the equation:

\[P_{t} =\frac{p^{t} \cdot q^{t}}{ \vec{b} \cdot q^{t}}\]

The price index values can be determined by normalizing the price levels by the first period as,

\[I_{t} = \frac{P_{t}}{P_{0}}\]

The caveat to using this matrix method is that the matrix $M = I_{N}-C+R$ should be invertible or equivalently $\det(M) \neq 0$. The matrix $M$ being singular is extremely rare and in such cases one may be tempted to use the Moore-Penrose pseudoinverse to compute the best fit solution for $\vec{b}$, but this will lead to a non-unique solution in some cases. Therefore the iterative method should be used as a fallback when the matrix $M$ is not invertible.

The code below is the implementation of the Geary-Khamis method using matrices. It is written in Python and uses the numpy and pandas library to solve the system of equations using the matrix method, but also the iterative method as a fallback. The matrix method is extremely efficient with a significantly shorter runtime.

from typing import List
import pandas as pd
import numpy as np

def diag(df: pd.DataFrame) -> pd.DataFrame:
    """
    Diagonal of a square dataframe, returned as a single-column frame so
    that it keeps its index for subsequent aligned division.
    """
    return pd.DataFrame(np.diag(df), index=df.index)

def geary_khamis(
    df: pd.DataFrame,
    price_col: str = 'price',
    quantity_col: str = 'quantity',
    date_col: str = 'month',
    product_id_col: str = 'id',
    method_type: str = 'matrix'
) -> List:
    """
    Geary-Khamis price index method

    Params:
        df: pd.DataFrame
            The dataframe containing the data.
        price_col: str
            The name of the price column.
        quantity_col: str
            The name of the quantity column.
        date_col: str
            The name of the date column.
        product_id_col: str
            The name of the product id column.
        method_type: str
            The type of method to use. Either 'matrix' or 'iterative'.

	Returns:
		List of price indices
	"""
    if method_type not in ('matrix', 'iterative'):
        raise ValueError('The method type must be `matrix` or `iterative`')

    # We need to deal with missing values and reshape the df for the
    # required vectors and matrices.
    df = df.dropna(how='all', axis=1).T.fillna(0)

    # Get number of unique products for the size of the vectors and
    # matrices.
    N = len(df.index.unique(level=product_id_col))

    # Matrices for the prices, quantities and weights.
    prices = df.loc[price_col]
    quantities = df.loc[quantity_col]
    weights = df.loc['weights']

    # Use iterative method directly if specified.
    if method_type == 'iterative':
        return _geary_khamis_iterative(prices, quantities)

    # Inverse of diagonal matrix with total quantities for each good over all
    # time periods as diagonal elements and matrix product of weights and
    # transpose of quantities to produce a square matrix, both required for C
    # matrix.
    q_matrix_inverse = np.diag(1/quantities.T.sum())
    prod_weights_qt_matrix = weights @ quantities.T

    # Product of above matrices to give the C square matrix, with the fixed
    # identity and R matrix, and c vector all required determine the quality
    # adjustment factors b.
    C_matrix = q_matrix_inverse @ prod_weights_qt_matrix
    R_matrix = np.zeros(shape=(N, N))
    R_matrix[:1] = 1
    
    # Define combo matrix used for isolating the singular matrices and
    # calculating the index values from the combo matrix `I_n - C + R`.
    combo_matrix = np.identity(N) - C_matrix + R_matrix
    
    if abs(np.linalg.det(combo_matrix)) <= 1e-7:
        # Fallback to iterative method for singular matrices.
        return _geary_khamis_iterative(prices, quantities)
    else:
        # Primary matrix method for non-singular matrices.
        return _geary_khamis_matrix(prices, quantities, combo_matrix)

def _geary_khamis_iterative(
    prices: pd.DataFrame,
    quantities: pd.DataFrame,
    no_of_iterations: int = 100,
    precision: float = 1e-8,
) -> pd.Series:
    """
    Geary-Khamis iterative method.
    
    Params:
        prices : pd.DataFrame
            The price dataframe.
        quantities : pd.DataFrame
            The quantity dataframe.
        no_of_iterations : int, defaults to 100
            The number of iterations to perform.
        precision : float, defaults to 1e-8
            The precision to use for the iterative method.

    Returns:
        pd.Series
            The price index values.
    """
    # Initialise index vals as 1's to find the solution with iteration.
    price_levels = pd.Series(1.0, index=prices.columns)
    quantity_share = quantities.T / quantities.sum(axis=1)

    # Iterate until we reach the set level of precision, or after a set
    # number of iterations if they do not converge.
    for _ in range(no_of_iterations):
        # Obtain matrices for iterative calculation.
        deflated_prices = prices / price_levels
        factors = diag(deflated_prices @ quantity_share)

        # Calculate new price levels from previous value.
        new_price_levels = (
            diag(prices.T @ quantities)
            .div(quantities.T @ factors)
            .squeeze()
        )

        pl_abs_diff = abs(price_levels - new_price_levels)
        
        if (pl_abs_diff <= precision).all():
            # Break loop when we reach given precision for final price levels.
            break
        else:
            # Otherwise set price level for next iteration.
            price_levels = new_price_levels

    # Normalize by first period for final output.
    return price_levels / price_levels.iloc[0]

def _geary_khamis_matrix(
    prices: pd.DataFrame,
    quantities: pd.DataFrame,
    combo_matrix: pd.DataFrame,
) -> pd.Series:
    """
    Geary-Khamis matrix method.
	
    Params:
        prices : pd.DataFrame
            The price dataframe.
        quantities : pd.DataFrame
            The quantity dataframe.
        combo_matrix : pd.DataFrame
            The combo matrix.
    
    Returns:
        pd.Series
            The price index values.
    """
    # Calculation of the vector b (factors) required to produce the
    # price levels. Corresponds to `b = [I_n - C + R]^-1 [1,0,..,0]^T`.
    # We use the Moore-Penrose inverse for the matrix inverse.
    factors = np.linalg.pinv(combo_matrix) @ np.eye(len(prices.index), 1)

    # Determine price levels to compute the final index values.
    price_levels = diag(prices.T @ quantities).div(quantities.T @ factors)

    # Normalize price levels to first period for final index values.
    index_vals = price_levels / price_levels.iloc[0]

    # Output as Pandas series for dynamic window.
    return index_vals.iloc[:, 0]

A worked implementation of this method, along with the other bilateral and multilateral methods (GEKS, Time Product Dummy, Time Dummy Hedonic, and the various bilateral index formulae), is available in my Python package PriceIndexCalc, which has a pandas backend and an experimental PySpark one:

pip install PriceIndexCalc

Which multilateral method won?

The question is no longer hypothetical, so it is worth recording where Geary-Khamis actually sits in 2026.

The ONS chose GEKS-Törnqvist on a 25-month window with a mean splice for groceries scanner data, and set out the reasoning in Introducing multilateral index methods into consumer price statistics. The short version:

  • Weighting philosophy. GEKS-Törnqvist is a superlative index: it uses the expenditure shares of both compared periods and makes no assumption about how consumers substitute between products. Geary-Khamis instead builds a single set of quality adjustment factors $b_n$ averaged across the whole window, which behaves like a quantity-weighted average and can be dominated by high-volume products.
  • Sensitivity to volume. Because $b_n$ is a quantity-weighted average price, GK is more exposed to products sold in very large volumes at low unit prices — a real problem in grocery data where pack sizes and units vary wildly.
  • International practice. GEKS-Törnqvist is what most early scanner-data adopters use (Australia, Canada, Belgium, Luxembourg, Norway, New Zealand), so the ONS choice keeps UK figures comparable.

None of that makes Geary-Khamis obsolete. It remains the standard method for spatial price comparisons — it is the backbone of the International Comparison Program and of purchasing power parity estimation, where the quality adjustment factors $b_n$ are exactly the “international prices” you want. It is also still used for temporal indices by some agencies, and it is a useful cross-check: if GK and GEKS-Törnqvist diverge sharply on the same window of data, that usually says something about product churn or unit-value effects in the underlying data rather than about the formulae.

The closed-form solution above is the reason GK stays practical at scale. On a 25-month window with tens of thousands of products, the iterative solution converges but you pay for it on every window in every product category, every month. Solving $\vec{b} = [I_N - C + R]^{-1}\vec{c}$ once is a single linear solve, and in the ONS pipeline work that difference was the gap between a method that could run within the production window and one that could not.

Two practical notes on the implementation, which have caught people out since:

  1. Do not silently substitute the Moore-Penrose pseudoinverse for the true inverse when $M$ is singular. It returns the minimum-norm solution, which is a solution but not necessarily the one satisfying the normalisation you imposed — hence the iterative fallback in the code above. (The _geary_khamis_matrix helper uses pinv for numerical robustness only after the determinant check has confirmed $M$ is well conditioned.)
  2. The determinant threshold of 1e-7 is scale-dependent. For very large $N$, prefer np.linalg.cond or an LinAlgError from np.linalg.solve as the singularity test, since $\det(M)$ underflows for large matrices regardless of conditioning.