Die Methode der kleinsten Quadrate ist ein Algorithmus aus dem Bereich des Machine Learnings.

Die Eingaben sind eine parametrisierte Funktion, die mehrerer Zahlen xi auf eine Zahl y = f(x1, ..., xn, p1, ... pn) abbildet, die Lernmenge aus vielen Tupeln aus zusammengehörigen xi und y und Startwerte für die Parameter pj.

Die Ausgabe besteht aus Parametern pj, für die die parametrisierte Funktion f anhand der xi den Wert y so gut wie möglich vorhersagt.

Der Güte der Vorhersage wird anhand der Summe der Quadrate der Differenzen von y und f(x1, ..., xn, p1, ... pn) bestimmt. Das übernimmt die Funktion residual mit den Hilfsfunktionen square und square-sum.

(defproc square (r) (* r r))

(defproc square-sum (rs s)
  (if
    (null? rs)
    s
    (square-sum (rest rs) (+ s (square (first rs))))))

(defproc residual (f ps xs ys)
  (square-sum
    (zip-with
      (lambda (x y) (- (f ps x) y))
      xs
      ys)
    0))

Das Ergebnis wird mit einem Gradientenverfahren berechnet. Der Gradient der Funktion f bezüglich der Parameter pj wird anhand von Sekantensteigungen angenähert. Dazu dienen die Funktionen secant, secants und gradient.

(defproc secant (f ph pt d xs ys)
  (/ (- (residual f (append ph (list (+ (first pt) d)) (rest pt)) xs ys)
        (residual f (append ph (list (- (first pt) d)) (rest pt)) xs ys))
     d
     2))

(defproc secants (f ph pt d xs ys)
  (if
    (null? pt)
    nil
    (cons
      (secant f ph pt d xs ys)
      (secants f (append ph (list (first pt))) (rest pt) d xs ys))))

(defproc gradient (f ps xs ys d)
  (secants f nil ps d xs ys))

Gradientenverfahren vermuten eine bessere Lösung entgegen der Richtung des Gradienten. Dabei ist die Schrittweite so zu bestimmen, dass sich tatsächlich eine Verbesserung ergibt.

(defproc improve (f ps xs ys gs r d b)
  (let
    ((hs
      (zip-with
        (lambda (p g) (- p (* b g)))
        ps
        gs)))
    (if
      (> r (residual f hs xs ys))
      (map-with (curry (approximate _ d)) hs)
      (improve f ps xs ys gs r d (/ b 3)))))

Die Funktion find-minimum führt die Verbesserungsschritte solange aus, bis eine vorgegebene Höchstzahl von Schritten erreicht wird oder der Gradient ein vorgegebenes Limit unterschreitet.

(defproc find-minimum (f ps xs ys d n)
  (let
    ((gs (gradient f ps xs ys (* d 1e-5))))
    (if
      (or (< (square-sum gs 0) (square d))
          (< n 1))
      ps
      (find-minimum
        f
        (improve f ps xs ys gs (residual f ps xs ys) d 1)
        xs
        ys
        d
        (- n 1)))))