12,507,674 members (51,799 online)

15.7K views
33 bookmarked
Posted

# Introduction to Numerical Solutions

, 24 Feb 2012 CPOL
An introduction to numerical solver algorithms with general purpose demonstration code.
 ```﻿/////////////////////////////////////////////////////////////////////////////// // // Solver.cs // // By Philip R. Braica (HoshiKata@aol.com, VeryMadSci@gmail.com) // // Distributed under the The Code Project Open License (CPOL) // http://www.codeproject.com/info/cpol10.aspx /////////////////////////////////////////////////////////////////////////////// // Using. using System; using System.Collections.Generic; using System.Linq; using System.Text; // Namespace namespace SolverDemo { /// /// Variate solver, supporting: /// Paris Genetic, /// Direct Search, /// Slope Variate. /// /// There are many other ways to solve things, steepest decent, /// Linear and non-linear least squares, weighted linear and non-linear /// least squares. Those are good approaches but other people have written /// about them, and explored them extensively, these are generally very useful, /// fast, and reliable alternatives I like. /// /// The Slope Variate solution is new, works to refine a solution, /// discovered by Phil Braica, 2011. /// public class Solver { /// /// Error computation delegate. /// For a given set of parameters, compare against known data, compute /// an error score where 0 = none, > 0 increasing in worseness. /// /// Parameters to try. /// Error public delegate double ComputeError_Delegate(double[] parameters); /// /// List of parameters to variate. /// public List Parameters = new List(); /// /// Step list. /// public List Steps = null; /// /// Best step to use. /// public Step BestStep { get; protected set; } /// /// Variate solver, paris style genetic algorithm. /// /// An N dimensional sampling of what the errors are over a region /// of N parameters is made, and the lowest is considered best. /// The best then seeds the next generation. /// /// It is up to the user to balance processing time vs. the likely smoothness /// of the error ... surface. For this reason the user selects the range /// for each parameter and not every variation is tested, the "center" /// maxStepsPerGeneration are used by sorting the possible things to try. /// /// The approach is more imune to settling on a local minimum instead of the /// best answer, but with each successive generation improves the answer /// more slowly than other approaches. Also if the user has set up the parameters /// just so, this approach may not converge to the ideal solution. /// /// Some solutions converge best if the software focuses on certain parameters first, /// then refines the others later. That "order of operations" quality is missing /// in this approach, which brings more "robustness" but causes less accuracy in /// some cases. /// /// Delegate to use. /// Number of generations /// Maximum number of test cases per generation. /// How much to expand search each generation, 1 = last step size, 1.5 = 1.5 * last step size. /// After this error amount, solution is good enough, return. /// The error from ComputeError for this solution. public double Solve_ParisGenetic( ComputeError_Delegate computeErrorFunction, int maxGenerations, int maxStepsPerGeneration, double retestFactor, double minStopError) { // No delegate? just return. if (computeErrorFunction == null) return -1; // Make a working list so we can adjust the parameters. List working = new List(); for (int i = 0; i < Parameters.Count; i++) { working.Add(new Parameter(Parameters[i])); } // Keep track of best score and best step. double bestScore = 0; Step bestStep = null; // For each generation. for (int i = 0; i < maxGenerations; i++) { // Make an initial step we will sweep through and variate // recursively. List steps = new List(); Step s = new Step(Parameters.Count); // Recursively compute the steps. // This is a similar algorithm to a multi-stride tree walk. computeSteps(steps, working, - 1, s); if (bestStep != null) { steps.Add(bestStep); // retry the best so new derivatives. } // Sort the steps so we only compute error terms for those // Closest to the center / likely solution if we need to // limit the number of steps we process. List sorted = new List(); sorted.AddRange(steps); sorted = sorted.OrderBy(x => x.CenterDistance).ToList(); // Compute errors for those we're interested, and keep track of what is best. bestStep = null; for (int j = 0; j < maxStepsPerGeneration && j < sorted.Count; j++) { sorted[j].Score = computeErrorFunction(sorted[j].ParamValues); sorted[j].Scored = true; if ((bestStep == null) || (sorted[j].Score < bestScore)) { bestStep = sorted[j]; bestScore = sorted[j].Score; } } Steps = steps; // If a best exists and it is good enough break. if ((bestStep != null) && (bestScore < minStopError)) { break; } // Adjust the step size for the next generation and run it. bool stepsChanged = false; for (int j = 0; j < working.Count; j++) { // Set new center. working[j].CenterValue = bestStep.ParamValues[j]; // Make sure initial is > 0 and > final, and if it was previously final, // no variation. double oldStep = working[j].InitialStepSize; double desiredStep = 0; if (oldStep != 0) { desiredStep = (working[j].MaxValue - working[j].MinValue) / (oldStep * oldStep); } if (working[j].InitialStepSize != working[j].FinalStepSize) { working[j].InitialStepSize = desiredStep; working[j].InitialStepSize = working[j].InitialStepSize < 0 ? 0 : working[j].InitialStepSize; working[j].InitialStepSize = working[j].InitialStepSize < working[j].FinalStepSize ? working[j].FinalStepSize : working[j].InitialStepSize; } stepsChanged |= (oldStep != working[j].InitialStepSize); // Adjust the bounds. double delta = working[j].InitialStepSize * retestFactor; delta = (oldStep != working[j].InitialStepSize) ? delta : 0; working[j].MinValue = working[j].CenterValue - delta; working[j].MaxValue = working[j].CenterValue + delta; } // If none of the steps changed then we're just retesting so break. if (stepsChanged == false) { break; } } // If we have a best set the solution. if (bestStep != null) { for (int i = 0; i < bestStep.ParamValues.Length; i++) { Parameters[i].Solution = bestStep.ParamValues[i]; } } BestStep = bestStep; // Return the score. return bestScore; } /// /// Variate solver, paris style genetic algorithm. /// /// The approach is the same as ParisGenetic, BUT multiple likely best cases /// are kept in each generation. /// /// Delegate to use. /// Number of generations /// Maximum number of test cases per generation. /// How much to expand search each generation, 1 = last step size, 1.5 = 1.5 * last step size. /// After this error amount, solution is good enough, return. /// The error from ComputeError for this solution. public double Solve_ParisGeneticMultiple( ComputeError_Delegate computeErrorFunction, int maxGenerations, int maxStepsPerGeneration, double retestFactor, double minStopError) { // No delegate? just return. if (computeErrorFunction == null) return -1; // Make a working list so we can adjust the parameters. List working = new List(); for (int i = 0; i < Parameters.Count; i++) { working.Add(new Parameter(Parameters[i])); } // Keep track of best score and best step. Step s = new Step(Parameters.Count); for (int i = 0; i < Parameters.Count; i++) { s.ParamValues[i] = Parameters[i].CenterValue; } List generation = new List(); generation.Add(s); int maxGenerationSize = 10; // For each generation. for (int i = 0; i < maxGenerations; i++) { // Make an initial step we will sweep through and variate // recursively. List nextGeneration = new List(); for (int j = 0; j < generation.Count; j++) { nextGeneration.AddRange(findBestGenetics( computeErrorFunction, working, generation[j], maxGenerationSize, maxStepsPerGeneration, retestFactor, minStopError)); } nextGeneration = nextGeneration.OrderBy(x => x.Score).ToList(); if (nextGeneration[0].Score < minStopError) { break; } if (nextGeneration.Count > maxGenerationSize) { nextGeneration.RemoveRange(maxGenerationSize, nextGeneration.Count - maxGenerationSize); } generation = nextGeneration; } for (int i = 0; i < generation[0].ParamValues.Length; i++) { Parameters[i].Solution = generation[0].ParamValues[i]; } return generation[0].Score; } protected List findBestGenetics( ComputeError_Delegate computeErrorFunction, List working, Step s, int maxGoodStepsPerGeneration, int maxStepsPerGeneration, double retestFactor, double minStopError) { List steps = new List(); computeStepsGenetic(steps, working, -1, s); // Sort the steps so we only compute error terms for those // Closest to the center / likely solution if we need to // limit the number of steps we process. List sorted = new List(); sorted.AddRange(steps); sorted = sorted.OrderBy(x => x.CenterDistance).ToList(); if (sorted.Count > maxStepsPerGeneration) { sorted.RemoveRange(maxStepsPerGeneration, sorted.Count - maxStepsPerGeneration); } for (int j = 0; j < sorted.Count; j++) { sorted[j].Score = computeErrorFunction(sorted[j].ParamValues); sorted[j].Scored = true; } sorted = sorted.OrderBy(x => x.Score).ToList(); if (sorted.Count > maxGoodStepsPerGeneration) { sorted.RemoveRange(maxGoodStepsPerGeneration, sorted.Count - maxGoodStepsPerGeneration); } return sorted; } /// /// Variate solver, Direct search aka. Hook-Jeves algorithm. /// /// Algorithm "walks" a parameter up and down a certain number of times attempting /// to reduce the error. Typcially each step divide the range of possibilities by 2. /// When this is done aggressively any ripples in the error vs parameter value might /// get skipped (either as a benefit or not) thus scaling the sub division and over /// testing the region is often good. /// /// Numerically this is faster than Direct Search, but it walks a single parameter /// at a time and thus not always converges, or converges instead to a local minimum. /// /// Delegate to use. /// Number of generations /// Maximum number of test cases per generation. /// How much to expand search each generation, 1 = cut step by half, generally between 0.75 and 1.5, above 2 and no effect. /// After this error amount, solution is good enough, return. /// The error from ComputeError for this solution. public double Solve_DirectSearch( ComputeError_Delegate computeErrorFunction, int maxGenerations, int maxStepsPerGeneration, double retestFactor, double minStopError) { // No delegate? just return. if (computeErrorFunction == null) return -1; retestFactor = retestFactor < 0.5 ? 0.5 : retestFactor > 1.999 ? 1.999 : retestFactor; // Make a working list we can adjust of the parameters. List working = new List(); for (int i = 0; i < Parameters.Count; i++) { working.Add(new Parameter(Parameters[i])); } // Keep track of best score and best step. Step bestStep = new Step(Parameters.Count); for (int i = 0; i < Parameters.Count; i++) { bestStep.ParamValues[i] = Parameters[i].CenterValue; } double bestScore = computeErrorFunction(bestStep.ParamValues); Step test = new Step(bestStep); retestFactor = retestFactor / 2; // For each generation. for (int i = 0; i < maxGenerations; i++) { // Select next parameter. for (int j = 0; j < Parameters.Count; j++) { for (int k = 0; k < maxStepsPerGeneration; k++) { Step lower = new Step(test); Step higher = new Step(test); lower.ParamValues[j] -= working[j].InitialStepSize; higher.ParamValues[j] += working[j].InitialStepSize; double lowerScore = computeErrorFunction(lower.ParamValues); double centerScore = computeErrorFunction(test.ParamValues); double higherScore = computeErrorFunction(higher.ParamValues); // If between lower and center, go to halfway between if ((lowerScore < centerScore) && (lowerScore < higherScore)) { test.ParamValues[j] -= working[j].InitialStepSize / 2; working[j].InitialStepSize *= retestFactor; } // Between higher and center go half between. if ((higherScore < centerScore) && (higherScore < lowerScore)) { test.ParamValues[j] += working[j].InitialStepSize / 2; working[j].InitialStepSize *= retestFactor; } // Between lower and upper, reduce step. if ((centerScore < lowerScore) && (centerScore < higherScore)) { working[j].InitialStepSize *= retestFactor; } if ((centerScore > lowerScore) && (centerScore > higherScore)) { if (lowerScore < higherScore) { test.ParamValues[j] -= working[j].InitialStepSize / 2; } else { test.ParamValues[j] += working[j].InitialStepSize / 2; } working[j].InitialStepSize *= retestFactor; } double err = computeErrorFunction(test.ParamValues); if (err < bestScore) { bestScore = err; bestStep = test; } if (bestScore < minStopError) { break; } } if (bestScore < minStopError) { break; } } if (bestScore < minStopError) { break; } } // If we have a best set the solution. if (bestStep != null) { for (int i = 0; i < bestStep.ParamValues.Length; i++) { Parameters[i].Solution = bestStep.ParamValues[i]; } } BestStep = bestStep; // Return the score. return bestScore; } /// /// Variate solver, slope relationship. /// /// This is simular to a direct search, but it makes certain assumptions. It is far less "robust", /// the initial guess needs to be relatively accurate, but it converges to an ideal solution faster than /// other approaches. /// /// At each stage: /// 1) Compute the curve that fits error vs. parameter k: err = F(k). /// 2) Compute best value of k. /// 3) Vary parameter j (using step 1 and 2 so that at each j value, the ideal k is used). /// 4) Compute the curve fit of jth parameter w/ ideal k @ j vs. error. /// 5) Compute the ideal jth parameter value. /// /// The output is an attempt to improve jth parameter by analyzing the way k interacts with j. /// The curve is always quadratic, so when the solution is close to being ideal, the lack of /// convergence is usually due to a shape that is either a cone or a sadle. Other solvers then spiral /// around the ideal solution but can't find it because they either under or overshoot parameter j, /// and then when refining parameter k have the same problem and it doesn't converge. /// This method variates slices through the curve to attempt to refine one of the parameters by /// actively refining the jth in terms of the kth parameter looking for a sadle point and driving to /// the minimum. /// /// /// Delegate to use. /// Number of generations /// Maximum number of test cases per generation. /// How much to expand search each generation, 1 = cut step by half, generally between 0.75 and 1.5, above 2 and no effect. /// After this error amount, solution is good enough, return. /// The error from ComputeError for this solution. public double Solve_SlopeVariate( ComputeError_Delegate computeErrorFunction, int maxGenerations, int maxStepsPerGeneration, double retestFactor, double minStopError) { // No delegate? just return. if (computeErrorFunction == null) return -1; retestFactor = retestFactor < 0.5 ? 0.5 : retestFactor > 1.999 ? 1.999 : retestFactor; // Make a working list we can adjust of the parameters. List working = new List(); for (int i = 0; i < Parameters.Count; i++) { working.Add(new Parameter(Parameters[i])); } // Keep track of best step. Step bestStep = new Step(Parameters.Count); for (int i = 0; i < Parameters.Count; i++) { bestStep.ParamValues[i] = Parameters[i].CenterValue; } bestStep.Score = computeErrorFunction(bestStep.ParamValues); bestStep.Scored = true; // Iterate through every pair of parameters maxGenerations number of times. for (int i = 0; i < maxGenerations; i++) { // Compute a nextBest. Step nextBest = new Step(bestStep); for (int j = 0; j < Parameters.Count; j++) { for (int k = 0; k < Parameters.Count; k++) { if (j == k) continue; // Test for j, k Step t = new Step(nextBest); // Vary j with k, the k with j. // By doing it twice, and not meeting a curve that doesn't fit, // two things occur: // 1) Instead of just refining j, and perhaps making the next // fit hard or impossible, j and k are addressed. // 2) We only consider worth continuing if there is a stability // of j on k and k on j. This ensures we didn't "improve" j // when really the curve fit was bad and the relationship // between j & k wasn't correct. for (int m = 0; m < maxStepsPerGeneration; m++) { // Improve jth of t, if no fit, don't improve t. t = vary_jk(computeErrorFunction, t, working, j, k); if (t == null) break; // Improve kth of t, t = vary_jk(computeErrorFunction, t, working, k, j); if (t == null) break; } if (t == null) continue; nextBest = t; nextBest.Score = computeErrorFunction(nextBest.ParamValues); nextBest.Scored = true; // If this variation is good enough, stop. if (nextBest.Score < minStopError) { break; } } // If this variation is good enough, stop. if (nextBest.Score < minStopError) { break; } } bestStep = nextBest; // If this variation is good enough, stop. if (bestStep.Score < minStopError) { break; } } BestStep = bestStep; BestStep.Scored = true; if (bestStep != null) { for (int i = 0; i < bestStep.ParamValues.Length; i++) { Parameters[i].Solution = bestStep.ParamValues[i]; } } // Return the score. return BestStep.Score; } /// /// Vary j, k, finding ideal kth params for each jth varied and then the ideal jth /// by minimal curve fit. This could be done recursively to compute an ideal solution. /// The problem is with measurement noise the likely stability as a general solution /// drops as the number of axises are increased. What is presented here is an iterative /// approach based on experiments. /// /// /// /// /// /// /// protected Step vary_jk(ComputeError_Delegate computeErrorFunction, Step s, List working, int j, int k) { int fp = 7; // fit points. double fp1 = fp - 1; double[] jx = new double[fp]; double[] jy = new double[fp]; double[] kx = new double[fp]; double[] ky = new double[fp]; Step vary = new Step(s); LinearLeastSquares lls = new LinearLeastSquares(); for (int ji = 0; ji < fp; ji++) { double fj = (fp1 - ((double)ji)) / fp1; double xfj = (working[j].MinValue * fj) + ((1 - fj) * working[j].MaxValue); jx[ji] = xfj; vary.ParamValues[j] = xfj; // Walk through kth at this j. for (int ki = 0; ki < fp; ki++) { double fk = (fp1 - ((double)ki)) / fp1; double xfk = (working[k].MinValue * fk) + ((1 - fk) * working[k].MaxValue); vary.ParamValues[k] = xfk; kx[ki] = xfk; ky[ki] = computeErrorFunction(vary.ParamValues); } // Now solve directly for what the zero point should be, given this j, varying k. lls.Solve_SecondOrder(kx, ky); // zy = a0 + a1*zx + a2*zx*zx // Slope is zero, is either a max or min, 2nd derivative // tells us if it's max or min. // zy' = a1 + 2*a2*zx = 0 -> a1 = -2*a2 *zx -> zx = -a1 / 2*a2 if (lls.Solution[2] <= 0) return null; double idealk = -0.5 * lls.Solution[1] / lls.Solution[2]; // idealk is the ideal value for the kth coefficient at this jth coefficient value. vary.ParamValues[k] = idealk; // This is now the smallest theoretical error when kth coefficient is perfectly // adjusted for jth coefficient. jy[ji] = computeErrorFunction(vary.ParamValues); } // Now the error equation of params j, k is jx, jy. // jx is walked through producing minimum errors for ideal k. lls.Solve_SecondOrder(jx, jy); // The minimum is at, slope = 0 if 2nd param is > 0. if (lls.Solution[2] <= 0) return null; Step t = new Step(s); double idealj = -0.5 * lls.Solution[1] / lls.Solution[2]; if (idealj < working[j].MinValue) return null; if (idealj > working[j].MaxValue) return null; t.ParamValues[j] = idealj; return t; } /// /// Compute the steps recursively starting at parameter p (-1 to init) based off /// of an initial step s (generated by this function when p = -1, then used recursively) /// based on the working parameters. /// /// /// /// /// protected void computeSteps(List steps, List working, int p, Step s) { // Done. Also acts as a "depth charge" for the recursion. if (p >= Parameters.Count) return; if (p == -1) { for (int i = 0; i < working.Count; i++) { s.ParamValues[i] = working[i].MinValue; } double sum = 0; for (int i = 0; i < working.Count; i++) { double px = s.ParamValues[i]; double dx = px - working[i].CenterValue; sum += dx * dx; } // This is an initial center distance, since we are sorting min to max // not computing or relying on a real "distance" from the center of the // solution space there is no need to compute the square root. s.CenterDistance = sum; steps.Add(s); p++; } // This call variates index p. double cx = working[p].CenterValue; Step s2 = null; for (double x = working[p].MinValue; x <= working[p].MaxValue; x += working[p].InitialStepSize) { s2 = new Step(s); s2.ParamValues[p] = x; double ox = s.ParamValues[p]; // Because no square root was done, we can subtract the old parameters distance // effect from the sum, and add in the new parameter values affect on the sum. s2.CenterDistance = s.CenterDistance + ((cx - x) * (cx - x)) - ((ox - cx) * (ox - cx)); steps.Add(s2); // Recurse. computeSteps(steps, working, p + 1, s2); } } /// /// Compute the steps recursively starting at parameter p (-1 to init) based off /// of an initial step s (generated by this function when p = -1, then used recursively) /// based on the working parameters. /// /// /// /// /// protected void computeStepsGenetic(List steps, List working, int p, Step s) { // s is the center. // Done. Also acts as a "depth charge" for the recursion. if (p >= Parameters.Count) return; if (p == -1) { // This is an initial center distance, since we are sorting min to max // not computing or relying on a real "distance" from the center of the // solution space there is no need to compute the square root. s.CenterDistance = 0; steps.Add(s); p++; } // This call variates index p. double cx = s.ParamValues[p]; Step s2 = null; double delta2 = 0.5 * (working[p].MaxValue - working[p].MinValue); for (double x = s.ParamValues[p] - delta2; x <= s.ParamValues[p] + delta2; x += working[p].InitialStepSize) { s2 = new Step(s); s2.ParamValues[p] = x; // Because no square root was done, we can subtract the old parameters distance // effect from the sum, and add in the new parameter values affect on the sum. s2.CenterDistance = s.CenterDistance + ((cx - x) * (cx - x)); steps.Add(s2); // Recurse. computeSteps(steps, working, p + 1, s2); } } /// /// Parameter class. /// public class Parameter { /// /// Copy constructor. /// /// public Parameter(Parameter p) { Name = p.Name; CenterValue = p.CenterValue; MaxValue = p.MaxValue; MinValue = p.MinValue; InitialStepSize = p.InitialStepSize; FinalStepSize = p.FinalStepSize; Solution = p.Solution; } /// /// Default constructor. /// public Parameter() { Name = ""; CenterValue = 0; MaxValue = 0; MinValue = 0; InitialStepSize = 0; FinalStepSize = 0; Solution = 0; } /// /// Constructor. /// /// /// /// /// /// /// /// public Parameter( string name, double centerValue, double maxValue, double minValue, double initialStepSize, double finalStepSize) { Name = name; CenterValue = centerValue; MaxValue = maxValue; MinValue = minValue; InitialStepSize = initialStepSize; Solution = 0; } /// /// Name. /// public string Name { get; set; } /// /// Center value (likely initial guess). /// public double CenterValue { get; set; } /// /// Max value. /// public double MaxValue { get; set; } /// /// Minimum value. /// public double MinValue { get; set; } /// /// Initial step size. /// public double InitialStepSize { get; set; } /// /// Final step size. /// public double FinalStepSize { get; set; } /// /// Final solution. /// public double Solution { get; set; } } /// /// A step. /// public class Step { /// /// Constructor. /// public Step() { } /// /// Constructor with size. /// /// public Step(int size) { ParamValues = new double[size]; } /// /// Copy constructor. /// /// public Step(Step s) { ParamValues = new double[s.ParamValues.Length]; for (int i = 0; i < s.ParamValues.Length; i++) { ParamValues[i] = s.ParamValues[i]; } Score = s.Score; } /// /// Scored flag. /// public bool Scored = false; /// /// Center distance. /// public double CenterDistance = 0; /// /// Parameter values. /// public double[] ParamValues; /// /// Score. /// public double Score; } } } ```

If a file you wish to view isn't highlighted, and is a text file (not binary), please let us know and we'll add colourisation support for it.

## Share

 Software Developer (Senior) KMC Systems United States
Phil is a Principal Software developer focusing on weird yet practical algorithms that run the gamut of embedded and desktop (PID loops, Kalman filters, FFTs, client-server SOAP bindings, ASIC design, communication protocols, game engines, robotics).

In his personal life he is a part time mad scientist, full time dad, and studies small circle jujitsu, plays guitar and piano.

## You may also be interested in...

 Pro Pro