DR error measures

This is a technical or mathematical page. While the subject may be of some relevance to music, the page treats the subject in technical language.

This article will describe several least-squares error measures for delta-rational chords. The idea behind least-squares error measures is to find the chord with the exact specified delta signature such that it deviates as little as possible from the chord that is meant to approximate that delta signature, and to measure the deviation. Least-squares error measures have the advantage of not fixing a particular interval in the chord when constructing the chord of best fit. However, like any other numerical measure of concordance or error, you should take them with a grain of salt.

Conventions and introduction

The idea motivating least-squares error measures on a chord as an approximation to a given delta signature is the following:

We want the error of a chord 1:r₁:r₂:...:r_n (in increasing order; we also take r₀ = 1), with n > 1, in the linear domain as an approximation to a delta-rational chord with signature +δ₁ +δ₂ ... +δ_n (possibly with some deltas free), i.e. the target chord is

$x:x+{\delta }_1:\cdots :x+\sum _{l=1}^n{\delta }_l$

with some root real-valued harmonic x. Let $D_0=0,D_i=\sum _{k=1}^i{\delta }_k$ be the delta signature +δ₁ +δ₂ ... +δ_n written cumulatively.

We want to measure the error without having to fix any dyad (as one might naively fix a dyad and measure errors in the other deltas in relation to the fixed dyad). To do this we solve a least-squares error problem: use a root-sum-square error function and optimize x (and any free deltas) to minimize that function.

Domain and error model

We have two choices:

• to measure either the linear (frequency ratio) error or the logarithmic (cents) one (called the domain);
• the collection of intervals to sum over (which we call the error model):
  • rooted: Only intervals from the root real-valued harmonic x are chosen.
  • pairwise: All intervals in the chord are compared.
  • all-steps: Only intervals between adjacent notes are compared.

We arrive at the following general formula: Let $[n]=\{0,1,2,...,n\},$ let $I\subseteq (\genfrac{}{}{0ex}{}{[n]}{2})$ be the error model, and let $f_D$ represent the domain function (identity for linear, or $\log$). Then the error function to be minimized by optimizing $x$ and any free deltas is:

$\sqrt{\sum _{i<j,\{i,j\}\in I}\left(f_D\right(\frac{x+D_j}{x+D_i})-f_D(\frac{r_j}{r_i}){)}^2}.$

Error function for various modes and error models

Domain	Error model	Error function
Linear	Rooted	$\sqrt{\sum _{i=1}^n(\frac{x+D_i}{x}-r_i{)}^2}=\sqrt{\sum _{i=1}^n(1+\frac{D_i}{x}-r_i{)}^2}$
	Pairwise	$\sqrt{\sum _{0\le i<j\le n}(\frac{x+D_j}{x+D_i}-\frac{r_j}{r_i}{)}^2}$
	All-steps	$\sqrt{\sum _{0\le i<n}(\frac{x+D_{i+1}}{x+D_i}-\frac{r_{i+1}}{r_i}{)}^2}$
Logarithmic (nepers)	Rooted	$\sqrt{\sum _{i=1}^n(\log \frac{x+D_i}{r_{i}x}{)}^2}$
	Pairwise	$\sqrt{\sum _{0\le i<j\le n}(\log \frac{x+D_j}{x+D_i}-\log \frac{r_j}{r_i}{)}^2}$
	All-steps	$\sqrt{\sum _{0\le i<n}(\log \frac{x+D_{i+1}}{x+D_i}-\log \frac{r_{i+1}}{r_i}{)}^2}$

To convert nepers to cents, multiply by $\frac{1200}{\log 2}.$

Solution methods

This section gets into the depths of mathematical optimization methods used to minimize DR error. (Optimization is a whole field and there are many different methods; the reason most of those methods exist is to find solutions when finding them symbolically is infeasible.)

Analytic

Fully DR, rooted linear

Setting the derivative to 0 gives us the closed-form solution

$x=\frac{\sum _{i=1}^n{D}_i}{-n+\sum _{i=1}^n{r}_i},$

which can be plugged back into

$\sqrt{\sum _{1=1}^n(1+\frac{D_i}{x}-r_i{)}^2}$

to obtain the least-squares linear error.

Partially DR (one related delta set, one free variable), rooted linear

Suppose we wish to approximate a target delta signature of the form $+{\delta }_1+?+{\delta }_3$ with the chord $1:r_1:r_2:r_3$ (where the +? is free to vary). The least-squares problem is to minimize by varying x and y the following:

$\sqrt{(\frac{x+{\delta }_1}{x}-r_1{)}^2+(\frac{x+{\delta }_1+y}{x}-r_2{)}^2+(\frac{x+{\delta }_1+y+{\delta }_3}{x}-r_3{)}^2},$

where y represents the free delta +?.

We can set the partial derivatives with respect to x and y of the inner expression equal to zero (since the derivative of sqrt() is never 0) and use SymPy to solve the system symbolically:

import sympy
x = sympy.Symbol("x", real=True)
y = sympy.Symbol("y", real=True)
d1 = sympy.Symbol("\\delta_{1}", real=True)
d2 = sympy.Symbol("\\delta_{2}", real=True)
d3 = sympy.Symbol("\\delta_{3}", real=True)
r1 = sympy.Symbol("r_1", real=True)
r2 = sympy.Symbol("r_2", real=True)
r3 = sympy.Symbol("r_3", real=True)
err_squared = ((x + d1) / x - r1) ** 2 + ((x + d1 + y) / x - r2) ** 2 + ((x + d1 + y + d3) / x - r3) ** 2
err_squared.expand()
err_squared_x = sympy.diff(err_squared, x)
err_squared_y = sympy.diff(err_squared, y)
sympy.nonlinsolve([err_squared_x, err_squared_y], [x, y])

The unique solution with x > 0 is

$x=\frac{2{\delta }_1+{\delta }_3+2y}{r_2+r_3-2},$

$y=\frac{-2{\delta }_1^2{r}_1+{\delta }_1^2{r}_2+{\delta }_1^2{r}_3-{\delta }_1{\delta }_3{r}_1+{\delta }_1{\delta }_3{r}_2-{\delta }_1{\delta }_3{r}_3+{\delta }_1{\delta }_3+{\delta }_3^2{r}_2-{\delta }_3^2}{2{\delta }_1{r}_1-2{\delta }_1-{\delta }_3{r}_2+{\delta }_3{r}_3}.$

Grid method (FDR case)

class GridSolution:
    def __init__(self, x, fx):
        self.x = x
        self.fx = fx
    
    def __repr__(self):
        return f"GridSolution(x={self.x:.5f}, fx={self.fx:.5f})"

def grid_method(f, window_size=100, coarse_steps=1000, fine_steps=1000):
    best_error = float("inf")
    best_x = 0
    
    coarse_step = window_size / coarse_steps
    # Fine step partitions one coarse step into smaller pieces
    fine_step = coarse_step / fine_steps 
    
    # --- Phase 1: Coarse Grid Search ---
    # Search in the window (0, window_size]
    for i in range(1, coarse_steps + 1):
        x = i * coarse_step
        fx = f(x)
        if fx < best_error:
            best_error = fx
            best_x = x

    # --- Phase 2: Fine Grid Search ---
    # Center the new search window around the best_x found in Phase 1
    # We search from (best_x - coarse_step/2) to (best_x + coarse_step/2)
    
    fine_window_lower = best_x - (coarse_step / 2)
    
    for j in range(1, fine_steps + 1):
        x = fine_window_lower + (j * fine_step)
        fx = f(x)
        if fx < best_error:
            best_error = fx
            best_x = x
            
    return GridSolution(best_x, best_error)

BFGS-B (one related delta set, arbitrary free deltas)

We let x₁ = x and include additional free variables x₂, ..., x_n, one for every additional +?, after coalescing segments of consecutive +?'s into one +? and after trimming leading and trailing free delta segments.

BFGS-B is a quasi-Newton optimization method (based on BFGS) particularly suited for this problem:

• The error function is smooth, allowing use of gradients
• Fast convergence, requiring at worst 20 iterations for accuracy
• Naturally deals with the x > 0 constraint using a log barrier and minimizing the transformed function using the unconstrained BFGS method
• Acceptable memory usage given a realistic number of parameters for practical DR chords (up to 3 interior free delta segments, thus 4 parameters).

It is a quasi-Newton method because it uses an approximation of the Hessian (matrix of mixed second partial derivatives) of the error function at each step.

In the Python implementation below, x represents the vector $(x_1,x_2,...,x_n),$ x0 is the initial guess for the solution, and f is the error function.

import numpy as np
import math

class BFGSSolution:
    def __init__(self, x, fx, iterations, success):
        self.x = x
        self.fx = fx
        self.iterations = iterations
        self.success = success

def numerical_gradient(f, x, eps=1e-8):
    grad = np.zeros_like(x)
    for i in range(len(x)):
        x_plus = x.copy()
        x_minus = x.copy()
        x_plus[i] += eps
        x_minus[i] -= eps
        grad[i] = (f(x_plus) - f(x_minus)) / (2 * eps)
    return grad

def bfgs(f, grad, x0, max_iterations=100, tolerance=1e-5):
    x = np.array(x0, dtype=float)
    n = len(x)
    fx = f(x)
    g = grad(x)
    
    # Approximate inverse Hessian
    # Initial approximation is the identity matrix
    H = np.eye(n)
    
    for i in range(max_iterations):
        # 0: Check convergence
        grad_norm = np.linalg.norm(g)
        if grad_norm < tolerance:
            return BFGSSolution(x, fx, i, True)
        
        # 1: Set search direction p (negative gradient direction transformed by H)
        # p = -H * g
        p = -H @ g
        
        # 2: Get alpha satisfying Wolfe conditions (Armijo rule and curvature condition)
        c1 = 1e-4
        c2 = 0.9
        max_line_search = 20
        rho_ls = 0.5
        
        gp = np.dot(g, p)
        alpha = 1.0
        
        for _ in range(max_line_search):
            x_next_guess = x + alpha * p
            
            # Check Armijo rule
            # f(x + alpha*p) <= f(x) + c1 * alpha * p^T * g
            if f(x_next_guess) <= fx + c1 * alpha * gp:
                # Check Curvature condition
                # -p^T * grad(x_next) <= -c2 * p^T * g
                g_next_guess = grad(x_next_guess)
                if -np.dot(p, g_next_guess) <= -c2 * gp:
                    break
            
            alpha *= rho_ls
        
        # 3: Set s = alpha * p and x_next = x + s
        s = alpha * p
        x_next = x + s
        
        # 4: Set y = grad(x_next) - grad(x)
        # We re-calculate grad(x_next) here to match the strict logic flow, 
        # though optimization could reuse g_next_guess from the successful line search.
        g_next = grad(x_next)
        y = g_next - g
        
        # 5: BFGS Update
        # Update H += U + V
        sy = np.dot(s, y)
        
        # Prevent division by zero if step size was extremely small
        if sy == 0: 
            # In a robust implementation, you might reset H to Identity here
            break 
            
        Hy = H @ y
        
        # Calculate scalar for the first term: (s^T y + y^T H y) / (s^T y)^2
        scalar1 = (sy + np.dot(y, Hy)) / (sy * sy)
        U = scalar1 * np.outer(s, s)
        
        # Calculate the second term matrices
        # W = (H y) s^T + s (y^T H)
        # Note: Since H is symmetric, y^T H is equivalent to (H y)^T
        W = np.outer(Hy, s) + np.outer(s, Hy)
        V = (-1 / sy) * W
        
        H = H + U + V
        
        # Update x, fx, and g for next iteration
        x = x_next
        fx = f(x_next)
        g = g_next

    return BFGSSolution(x, fx, max_iterations, False)

def bfgs_barrier(f, bounds, x0, history_size=10, max_iterations=100, tolerance=1e-5, barrier_weight=1e-4):
    """
    Solves bounded optimization using a Log-Barrier method wrapped around BFGS.
    Note: x0 must be strictly feasible (inside bounds).
    """
    def transformed_f(x):
        penalty = 0
        for i, val in enumerate(x):
            lower, upper = bounds[i]
            
            # Hard cutoff prevents log domain errors during line search exploration
            if (lower is not None and val <= lower) or (upper is not None and val >= upper):
                return float("inf")
            
            # Log barrier penalties
            if lower is not None:
                penalty -= barrier_weight * math.log(val - lower)
            if upper is not None:
                penalty -= barrier_weight * math.log(upper - val)
        
        return f(x) + penalty
    
    # Use the transformed function for gradients as well
    grad = lambda x: numerical_gradient(transformed_f, x)
    
    result = bfgs(transformed_f, grad, x0, history_size, max_iterations, tolerance)
    
    # Restore actual function value
    result.fx = f(result.x)
    return result

L-BFGS-B (one related delta set, arbitrary free deltas)

L-BFGS-B is an approximation to BFGS-B limiting memory usage, particularly suited for high-dimensional problems but nevertheless agreeing very well with BFGS-B for typical DR test cases. Python code for the L-BFGS method is provided below.

import numpy as np
import math

class LBFGSSolution:
    def __init__(self, x, fx, iterations, success):
        self.x = x
        self.fx = fx
        self.iterations = iterations
        self.success = success

def numerical_gradient(f, x, eps=1e-8):
    grad = np.zeros_like(x)
    for i in range(len(x)):
        x_plus = x.copy()
        x_minus = x.copy()
        x_plus[i] += eps
        x_minus[i] -= eps
        grad[i] = (f(x_plus) - f(x_minus)) / (2 * eps)
    return grad

def l_bfgs(f, grad, x0, history_size=10, max_iterations=100, tolerance=1e-5):
    x = np.array(x0, dtype=float) # Ensure float array
    fx = f(x)
    g = grad(x)
    
    s_history = []
    y_history = []
    rho_history = []
    
    for iteration in range(max_iterations):
        grad_norm = np.linalg.norm(g)
        if grad_norm < tolerance:
            return LBFGSSolution(x, fx, iteration, True)
        
        # --- Two-loop recursion ---
        q = g.copy()
        
        # We need to store alpha to use it in the second loop
        m = len(s_history)
        alpha = [0.0] * m 
        
        # First loop (Backward)
        for i in range(m - 1, -1, -1):
            alpha[i] = rho_history[i] * np.dot(s_history[i], q)
            q = q - alpha[i] * y_history[i]
        
        # Initial Hessian approximation scaling
        gamma = 1.0
        if m > 0:
            gamma = np.dot(s_history[-1], y_history[-1]) / np.dot(y_history[-1], y_history[-1])
        
        z = q * gamma

        # Second loop (Forward)
        for i in range(m):
            beta = rho_history[i] * np.dot(y_history[i], z)
            z = z + s_history[i] * (alpha[i] - beta)
        
        p = -z # Search direction

        # --- Line search (Backtracking) ---
        step_size = 1.0
        c1 = 1e-4
        rho_ls = 0.5
        max_line_search = 20
        
        # Check for sufficient decrease (Armijo rule)
        # f(x + alpha*p) <= f(x) + c1 * alpha * p^T * g
        
        g_dot_p = np.dot(g, p)
        
        for _ in range(max_line_search):
            x_new = x + step_size * p
            fx_new = f(x_new)
            
            if fx_new <= fx + c1 * step_size * g_dot_p:
                break
            step_size *= rho_ls
        else:
            # If line search fails to find a better point, stop or accept best attempt
            # (Here we just accept the last reduced step_size)
            pass 
        
        # --- Update History ---
        s = x_new - x
        g_new = grad(x_new)
        y = g_new - g
        
        sy = np.dot(s, y)
        
        if sy > 1e-10: # Ensure curvature condition
            s_history.append(s)
            y_history.append(y)
            rho_history.append(1.0 / sy)
            
            # Maintain history size
            if len(s_history) > history_size:
                s_history.pop(0)
                y_history.pop(0)
                rho_history.pop(0)
            
        x = x_new
        fx = fx_new
        g = g_new
    
    return LBFGSSolution(x, fx, max_iterations, False)

External links

• Inthar's DR chord explorer (includes least-squares error calculation in both domains and multiple error models)