Remez Algorithm - Source Code
This page displays the source code of key functions and methods used in the Remez algorithm implementation.
initializeCFRemezNodes(targetFunc, a, b, numNodes) {
if (!Number.isFinite(a) || !Number.isFinite(b) || a >= b || numNodes < 2) {
throw new Error(`Invalid inputs: a=${a}, b=${b}, numNodes=${numNodes}`);
}
const n = numNodes - 1; // degree of approximation
const numSamples = Math.max(10 * n, 127); // More samples for stability
console.log('numSamples:', numSamples, 'numNodes:', numNodes, 'degree:', n);
// Step 1: Compute Chebyshev coefficients more robustly
const chebCoeffs = new Array(2 * n + 1).fill(0);
// Use Chebyshev-Gauss-Lobatto nodes for better accuracy
for (let k = 0; k < numSamples; k++) {
const t = Math.cos(Math.PI * k / (numSamples - 1)); // t in [-1, 1]
const x = (b - a) * (t + 1) / 2 + a; // map to [a, b]
const y = targetFunc(x);
if (!Number.isFinite(y)) {
console.warn(`Non-finite function value at x=${x}`);
continue;
}
// Compute Chebyshev coefficients using discrete cosine transform
for (let j = 0; j < 2 * n + 1 && j < numSamples; j++) {
const Tj = this.chebyshevPolynomial(j, t);
chebCoeffs[j] += y * Tj;
}
}
// Proper DCT normalization for Chebyshev coefficients
for (let j = 0; j < chebCoeffs.length; j++) {
if (j === 0) {
chebCoeffs[j] /= numSamples; // c_0 normalization
} else {
chebCoeffs[j] *= 2.0 / numSamples; // c_k normalization for k > 0
}
// Clean up small coefficients
if (Math.abs(chebCoeffs[j]) < 1e-14) {
chebCoeffs[j] = 0;
}
}
console.log('Chebyshev coefficients:', chebCoeffs.slice(0, Math.min(10, chebCoeffs.length)));
// Step 2: Set up CF approximation problem
// For CF approximation, we need to solve a different eigenvalue problem
// We construct a Toeplitz matrix from the Chebyshev coefficients
let cfCoefficients;
let cfApprox;
try {
// Method 1: Direct coefficient truncation (simpler, often more stable)
cfCoefficients = chebCoeffs.slice(0, n + 1);
// Create the CF approximation function
cfApprox = (x) => {
// Transform x from [a,b] to [-1,1]
const t = 2 * (x - a) / (b - a) - 1;
let result = cfCoefficients[0]; // T_0 = 1, use coefficient directly
for (let k = 1; k < cfCoefficients.length; k++) {
result += cfCoefficients[k] * this.chebyshevPolynomial(k, t);
}
return result;
};
console.log('CF coefficients (truncated):', cfCoefficients);
} catch (error) {
console.error('CF coefficient computation failed:', error);
// Fallback: Use Chebyshev interpolation nodes
const chebNodesFallback = Array.from(
{ length: numNodes },
(_, k) => (b - a) * (Math.cos(Math.PI * k / (numNodes - 1)) + 1) / 2 + a
).sort((x, y) => x - y);
return {
coeffs: new Array(n + 1).fill(0),
nodes: chebNodesFallback,
cfApprox: x => targetFunc(x),
maxError: Infinity
};
}
// Step 3: Find alternating extrema more robustly
const errorFunc = (x) => {
const target = targetFunc(x);
const approx = cfApprox(x);
return target - approx;
};
// Sample the error function to check quality
/*const testPoints = remMath.linspace(a, b, 50);
const errors = testPoints.map(x => ({
x: x,
target: targetFunc(x),
approx: cfApprox(x),
error: errorFunc(x)
}));
console.log('Error function samples (first 5):', errors.slice(0, 5));
console.log('Max error in sample:', Math.max(...errors.map(e => Math.abs(e.error))));*/
// Find extrema with better initial guess
let extrema;
try {
extrema = this.findAlternatingExtrema(targetFunc, cfApprox, a, b);
console.log('Found extrema:', extrema);
} catch (error) {
console.error('Extrema finding failed:', error);
// Use evenly spaced points as fallback
extrema = Array.from({ length: numNodes }, (_, k) => {
const x = a + k * (b - a) / (numNodes - 1);
return { x: x, y: errorFunc(x) };
});
}
const nodes = extrema.map(e => e.x);
const maxError = Math.max(...extrema.map(e => Math.abs(e.y)));
return {
coeffs: cfCoefficients,
nodes: nodes,
cfApprox: cfApprox,
maxError: maxError
};
}initializeL2RemezNodes(targetFunc, a, b, numNodes) {
const l2A = remMath.polynomialL2Approximation(targetFunc, numNodes-2, a, b);
const lsApprox = l2A.lsApprox;
let extrema = this.findAlternatingExtrema(targetFunc, l2A.lsApprox, a, b, numNodes);
const nodes = extrema.map(function(e) { return e.x; });
return {
coeffs: l2A.coeffs,
nodes: nodes,
lsApprox: l2A.lsApprox
};
}initializeDiscreteRemezNodes(targetFunc, a, b, numNodes, gridDensity = 100) {
const degree = numNodes - 2;
const numGridPoints = Math.max(1000, gridDensity * numNodes);
const gridPoints = Array.from({ length: numGridPoints }, (_, i) =>
a + i * (b - a) / (numGridPoints - 1)
);
const targetValues = gridPoints.map(x => targetFunc(x));
const discreteResult = this.solveDiscreteRemez(gridPoints, targetValues, degree);
const extrema = this.findDiscreteExtremaFromSolution(
gridPoints,
targetValues,
discreteResult.approxValues,
numNodes
);
const nodes = extrema.map(e => e.x);
return {
coeffs: discreteResult.coeffs,
nodes: nodes,
discreteApprox: discreteResult.approxFunc,
maxError: Math.max(...extrema.map(e => Math.abs(e.error)))
};
}findReferenceSimple(errorFunction, coefficients, currentReference) {
const coarseSteps = parseInt(document.getElementById('bruteSteps').value);
const refinementSteps = 10;
const refinementIterations = 5;
const eps = 0.000001/n;
let newReference = [];
for (let k = 0; k < currentReference.length; k++) {
const left = (k==0) ? aP.intervalStart : currentReference[k-1].x + eps;
const right = (k==currentReference.length-1) ? aP.intervalEnd : currentReference[k+1].x;
// Safety check: ensure valid interval
if (left >= right || !isFinite(left) || !isFinite(right)) {
console.warn(`Invalid interval at k=${k}: [${left}, ${right}]`);
newReference.push({...currentReference[k], originalIndex: k});
continue;
}
// Initial coarse search
let extremum = currentReference[k];
let delta = (right - left) / coarseSteps;
// Safety check: ensure valid delta
if (delta <= 0 || !isFinite(delta)) {
console.warn(`Invalid delta at k=${k}: ${delta}`);
newReference.push({...extremum, originalIndex: k});
continue;
}
// Use step counter instead of relying only on currentX <= right
for (let step = 0; step <= coarseSteps; step++) {
const currentX = left + step * delta;
if (currentX > right) break; // Additional safety
const currentY = errorFunction(coefficients, currentX);
if (!isFinite(currentY)) continue; // Skip invalid results
if ((currentReference[k].y > 0 && currentY > extremum.y) ||
(currentReference[k].y < 0 && currentY < extremum.y)) {
extremum = {x: currentX, y: currentY};
}
}
// Refinement phase
for (let i = 0; i < refinementIterations; i++) {
const refinementLeft = Math.max(left, extremum.x - delta);
const refinementRight = Math.min(right, extremum.x + delta);
delta = (refinementRight - refinementLeft) / refinementSteps;
// Safety check for refinement delta
if (delta <= 0 || !isFinite(delta)) break;
// Use step counter for refinement too
for (let step = 0; step <= refinementSteps; step++) {
const currentX = refinementLeft + step * delta;
if (currentX > refinementRight) break;
const currentY = errorFunction(coefficients, currentX);
if (!isFinite(currentY)) continue;
if ((currentReference[k].y > 0 && currentY > extremum.y) ||
(currentReference[k].y < 0 && currentY < extremum.y)) {
extremum = {x: currentX, y: currentY};
}
}
}
newReference.push({...extremum, originalIndex: k});
}
return newReference;
}findReferenceParabolic(calculateError, optimizationCoefficients, currentReferencePoints) {
const maxInterpolationIterations = parseInt(document.getElementById('parabolicSteps').value);
const minimumPointSeparation = 0.01;
const convergenceThreshold = 1e-10;
const epsilon = 1e-12;
const goldenRatio = 0.618033988749;
// Enhanced parabolic interpolation with better numerical stability
function interpolateParabola(x0, x1, x2, y0, y1, y2) {
// Use more numerically stable form
const h1 = x1 - x0;
const h2 = x2 - x1;
const d1 = (y1 - y0) / h1;
const d2 = (y2 - y1) / h2;
const denominator = h1 + h2;
if (Math.abs(denominator) < epsilon || Math.abs(d2 - d1) < epsilon) {
return (x0 + x2) / 2; // Fallback to midpoint
}
const a = (d2 - d1) / denominator;
if (Math.abs(a) < epsilon) {
return (x0 + x2) / 2; // Nearly linear, use midpoint
}
// Compute vertex of parabola
const vertex = x1 - d1 / (2 * a) - h1 / 2;
return vertex;
}
// Adaptive bracket width calculation
function calculateInitialBracket(pointIndex, currentPoints) {
const baseWidth = 0.1; // Default bracket width
let adaptiveWidth = baseWidth;
if (pointIndex > 0) {
const leftSpacing = currentPoints[pointIndex].x - currentPoints[pointIndex - 1].x;
adaptiveWidth = Math.min(adaptiveWidth, leftSpacing * 0.3);
}
if (pointIndex < currentPoints.length - 1) {
const rightSpacing = currentPoints[pointIndex + 1].x - currentPoints[pointIndex].x;
adaptiveWidth = Math.min(adaptiveWidth, rightSpacing * 0.3);
}
return Math.max(adaptiveWidth, minimumPointSeparation * 2);
}
// Enhanced caching with LRU eviction
class ErrorCache {
constructor(maxSize = 1000) {
this.cache = new Map();
this.maxSize = maxSize;
}
get(x) {
const key = x.toFixed(12); // Round to avoid floating point issues
if (this.cache.has(key)) {
const value = this.cache.get(key);
this.cache.delete(key);
this.cache.set(key, value); // Move to end (most recent)
return value;
}
const value = calculateError(optimizationCoefficients, x);
if (this.cache.size >= this.maxSize) {
const firstKey = this.cache.keys().next().value;
this.cache.delete(firstKey);
}
this.cache.set(key, value);
return value;
}
}
const errorCache = new ErrorCache();
// Main optimization loop
let optimizedReferencePoints = [];
for (let pointIndex = 0; pointIndex < currentReferencePoints.length; pointIndex++) {
const currentPoint = currentReferencePoints[pointIndex];
// Calculate boundaries
const minAllowedX = (pointIndex === 0) ?
aP.intervalStart :
optimizedReferencePoints[pointIndex - 1].x + minimumPointSeparation;
const maxAllowedX = (pointIndex === currentReferencePoints.length - 1) ?
aP.intervalEnd :
currentReferencePoints[pointIndex + 1].x - minimumPointSeparation;
// Initialize with adaptive bracket
const bracketWidth = calculateInitialBracket(pointIndex, currentReferencePoints);
let leftX = Math.max(minAllowedX, currentPoint.x - bracketWidth);
let midX = currentPoint.x;
let rightX = Math.min(maxAllowedX, currentPoint.x + bracketWidth);
// Ensure proper ordering and spacing
if (rightX - leftX < 2 * minimumPointSeparation) {
const center = (leftX + rightX) / 2;
leftX = Math.max(minAllowedX, center - minimumPointSeparation);
rightX = Math.min(maxAllowedX, center + minimumPointSeparation);
midX = center;
}
let bestX = midX;
let bestY = errorCache.get(midX);
let iterationCount = 0;
let consecutiveSmallSteps = 0;
let lastStepSize = Infinity;
// Target direction: maximize absolute error
const targetSign = Math.sign(currentPoint.y);
while (iterationCount < maxInterpolationIterations) {
const leftY = errorCache.get(leftX);
const midY = errorCache.get(midX);
const rightY = errorCache.get(rightX);
// Try parabolic interpolation
let newX = interpolateParabola(leftX, midX, rightX,
leftY * targetSign, midY * targetSign, rightY * targetSign);
// Clamp to valid range
newX = Math.max(minAllowedX, Math.min(maxAllowedX, newX));
// Avoid getting stuck at boundaries
if (newX === minAllowedX || newX === maxAllowedX) {
newX = (leftX + rightX) / 2;
}
const stepSize = Math.abs(newX - bestX);
// Check convergence
if (stepSize < convergenceThreshold) {
consecutiveSmallSteps++;
if (consecutiveSmallSteps >= 3) break;
} else {
consecutiveSmallSteps = 0;
}
// Detect oscillation or stagnation
if (stepSize > lastStepSize * 0.9 && iterationCount > 3) {
// Switch to golden section if parabolic isn't converging
const range = rightX - leftX;
if (newX > midX) {
newX = midX + goldenRatio * (rightX - midX);
} else {
newX = midX - goldenRatio * (midX - leftX);
}
newX = Math.max(minAllowedX, Math.min(maxAllowedX, newX));
}
const newY = errorCache.get(newX);
// Update best point if we found improvement
if (Math.abs(newY) > Math.abs(bestY)) {
bestX = newX;
bestY = newY;
}
// Update bracket based on function behavior
if (newX < midX) {
if (Math.abs(newY) > Math.abs(midY)) {
rightX = midX;
midX = newX;
} else {
leftX = newX;
}
} else if (newX > midX) {
if (Math.abs(newY) > Math.abs(midY)) {
leftX = midX;
midX = newX;
} else {
rightX = newX;
}
}
// Ensure bracket doesn't collapse
if (rightX - leftX < 2 * convergenceThreshold) {
const expansion = Math.max(convergenceThreshold * 10, (maxAllowedX - minAllowedX) * 0.01);
const center = (leftX + rightX) / 2;
leftX = Math.max(minAllowedX, center - expansion);
rightX = Math.min(maxAllowedX, center + expansion);
midX = center;
}
lastStepSize = stepSize;
iterationCount++;
}
// Final adjustment for minimum separation
if (optimizedReferencePoints.length > 0) {
const minRequiredX = optimizedReferencePoints[optimizedReferencePoints.length - 1].x + minimumPointSeparation;
if (bestX < minRequiredX) {
bestX = minRequiredX;
bestY = errorCache.get(bestX);
}
}
optimizedReferencePoints.push({
x: bestX,
y: bestY,
iterations: iterationCount,
converged: consecutiveSmallSteps >= 3 || iterationCount < maxInterpolationIterations
});
}
console.log("Optimized Reference Points:", optimizedReferencePoints);
return optimizedReferencePoints;
}