Discrete best uniform approximation with Remez algorithm

remez.py

import numpy as np
import scipy.signal

def discrete_remez(f, V, tol=1e-12, maxiter=1000, startnodes=None):
    """Best uniform approximation of the vector f by a linear combination of the columns of V.
    Uses a discrete Remez algorithm with single point exchange.
    
    If it converges, this function returns a vector c such that V.dot(c) is the best uniform
    approximation to f.
    """
    # orthonormalize
    V, Lam, W = np.linalg.svd(V, full_matrices=False)
    k = V.shape[1]   # dimension of space
    
    if startnodes is not None:
        ii = np.array(startnodes).ravel()
        if len(ii) != k + 1:
            raise ValueError('invalid number of starting nodes: %d, should be %d' % (len(ii), k+1))
    else:
        # find initial approximation by least squares
        c = V.T.dot(f) #np.linalg.solve(V.T.dot(V), V.T.dot(f)) # non-ortho case
        # try to find peaks of error function
        e = abs(f - V.dot(c))
        e0 = np.pad(e, [1,1], mode='constant')  # pad with 0 to find boundary maxima
        ii = scipy.signal.argrelmax(e0)[0] - 1  # subtract padded index
        if len(ii) != k + 1:
            raise RuntimeError('wrong number of initial peaks found: %d' % len(ii))
            
    s1 = (-1) ** np.arange(k+1)
    
    for it in range(maxiter):
        # solve leveling equation
        sol = np.linalg.solve(np.hstack((V[ii], s1[:,None])), f[ii])
        c, lam = sol[:-1], sol[-1]
        
        # check error - does it equioscillate?
        e = abs(f - V.dot(c))
        delta = max(e) - min(e[ii])
        if delta < tol:
            return W.T.dot(c / Lam)
        
        # single point exchange
        imax = np.argmax(e)
        # find closest node
        j = np.argmin(abs(imax - ii))
        # replace it
        ii[j] = imax

    print('did not converge: delta = %e' % delta)
    return W.T.dot(c / Lam)