Parameteridentifikation

Motivation

Den Anstoß für diese Erweiterung gab eine praktische Frage aus dem wissenschaftlichen Rechnen: Wie lässt sich ein vorhandenes numerisches Modell zwischen MATLAB, Python und Fortran übertragen, ohne seine mathematische Aussage oder seine numerische Qualität zu verändern?

Eine solche Portierung ist keine reine Sprachfrage. Zuerst muss geklärt werden, ob beide Programme wirklich dasselbe mathematische Modell, dieselbe Diskretisierung und dieselben Rand- beziehungsweise Quadraturbehandlungen verwenden. Erst danach sind Aussagen über numerische Übereinstimmung, Wartbarkeit oder Geschwindigkeit sinnvoll.

Das Modellproblem dieses Beitrags bietet dafür einen überschaubaren Referenzfall. Ausgehend vom vorangehenden Sage-Skript diskretisieren die MATLAB-, Python- und Fortran-Fassungen denselben Integraloperator und lösen dieselben direkten beziehungsweise regularisierten Gleichungssysteme. Der Vergleich konzentriert sich auf zwei Fragen: Stimmen die numerischen Ergebnisse überein, und wie unterscheiden sich Laufzeit und Speicherbedarf der konkreten Implementierungen?

Für kleine und mittlere Diskretisierungen steht die mathematische Validierung im Vordergrund. Große Matrizen werden ausschließlich als Performance- und Skalierungstest verwendet, da die direkte inverse Aufgabe mit wachsendem n extrem schlecht konditioniert ist.

Gemeinsame diskrete Referenz

Alle drei Implementierungen verwenden die äquidistanten Stützstellen \(x_i=i/n\), \(i=0,\ldots,n\), und die Gewichte der summierten Trapezregel. Mit

\[
K_{ij}=\exp(x_i x_j),
\qquad
W=\operatorname{diag}(w_0,\ldots,w_n),
\qquad
A=KW
\]

lautet die direkte diskrete Aufgabe

\[
A\boldsymbol{\alpha}=\boldsymbol{y}.
\]

Die konsistent diskretisierte Tikhonov-Gleichung wird in allen Programmen in derselben Form gelöst:

\[
(\beta I+A A)\boldsymbol{\alpha}_{\beta}=A\boldsymbol{y}.
\]

Damit vergleichen wir tatsächlich dieselbe mathematische Aufgabe. Insbesondere darf die diskrete adjungierte Abbildung hier nicht stillschweigend durch die gewöhnliche transponierte Matrix ersetzt werden; die Quadraturgewichte müssen konsistent berücksichtigt bleiben.

MATLAB-Implementierung

Die MATLAB-Fassung bildet Stützstellen, Quadraturgewichte und Matrix vektorisiert und verwendet die lineare Algebra des installierten MATLAB-Backends. Direkte und regularisierte Lösung werden getrennt ausgewertet, damit Konditionierung, Rekonstruktionsfehler und Laufzeit nicht miteinander vermischt werden.

% Diskretisierung und konsistente Tikhonov-Regularisierung
n = 32;
beta = 1.0e-6;

x = (0:n).'/n;
weights = ones(n+1, 1)/n;
weights([1, end]) = 1/(2*n);

A = exp(x*x.').*weights.';
y = expm1(x+1)./(x+1);
lambdaExact = exp(x);

conditionA = cond(A, 2);
alphaDirect = A\y;
directError = max(abs(alphaDirect-lambdaExact));

lhs = beta*eye(n+1) + A*A;
rhs = A*y;
alphaRegularized = lhs\rhs;
regularizedError = max(abs(alphaRegularized-lambdaExact));

Beim Matrixaufbau bezeichnet .* die elementweise Multiplikation des Kerns mit den Quadraturgewichten. Die Produkte A*A und A*y sind dagegen gewöhnliche Matrixprodukte.

Der vollständige MATLAB-Quellcode mit Tests und Benchmarkprogrammen ist auf GitHub verfügbar.

Python-/NumPy-Implementierung

Die Python-Fassung verwendet NumPy und folgt derselben Reihenfolge von Matrixaufbau, direkter Lösung und Tikhonov-Regularisierung. Für einen fairen Vergleich werden keine sprachspezifischen Vereinfachungen verwendet, die das mathematische Verfahren verändern würden.

# Diskretisierung und konsistente Tikhonov-Regularisierung
import numpy as np

n = 32
beta = 1.0e-6

x = np.arange(n + 1, dtype=np.float64) / n
weights = np.full(n + 1, 1.0 / n, dtype=np.float64)
weights[[0, -1]] = 1.0 / (2.0 * n)

A = np.exp(np.outer(x, x)) * weights[np.newaxis, :]
y = np.expm1(x + 1.0) / (x + 1.0)
lambda_exact = np.exp(x)

condition_A = np.linalg.cond(A)
alpha_direct = np.linalg.solve(A, y)
direct_error = np.max(np.abs(alpha_direct - lambda_exact))

lhs = beta * np.eye(n + 1) + A @ A
rhs = A @ y
alpha_regularized = np.linalg.solve(lhs, rhs)
regularized_error = np.max(np.abs(alpha_regularized - lambda_exact))

NumPy verwendet * für die elementweise Multiplikation und @ für Matrixprodukte. Abgesehen von diesen syntaktischen Unterschieden führen beide Listings dieselben Rechenschritte in derselben Reihenfolge aus.

Der vollständige Python-/NumPy-Quellcode mit Tests und Benchmarkprogrammen ist auf GitHub verfügbar.

Fortran-90/95-Implementierung

Die Fortran-Fassung ist eine eigenständige Referenzimplementierung desselben Modells. Der numerische Kern ist in Module für Präzision, Modellaufbau, lineare Algebra und Auswertung gegliedert und wird mit strengem Fortran 95 übersetzt. DGESV, DGESVD, DGEMM und DGEMV binden die etablierten BLAS-/LAPACK-Routinen direkt ein. Nur das getrennte Benchmark-Modul nutzt für die hochauflösende Wandzeitmessung eine dokumentierte neuere SYSTEM_CLOCK-Erweiterung von gfortran.

Der folgende Ausschnitt stammt aus dem tatsächlichen Modellmodul. Er zeigt die konkrete Datenstruktur, ihre explizite Speicherverwaltung und den Aufbau der diskreten Aufgabe mit Fortran-95-Mitteln:

module gwr_model
    use gwr_precision, only: dp
    implicit none

    type, public :: problem_data
        integer :: n = 0
        real(dp), pointer :: x(:)            => null()
        real(dp), pointer :: weights(:)      => null()
        real(dp), pointer :: matrix(:, :)    => null()
        real(dp), pointer :: data(:)         => null()
        real(dp), pointer :: exact_lambda(:) => null()
    end type problem_data

contains

    subroutine build_problem(n, problem, info)
        integer, intent(in) :: n
        type(problem_data), intent(out) :: problem
        integer, intent(out) :: info
        integer :: i, j, stat

        info = 0
        nullify(problem%x, problem%weights, problem%matrix, &
                problem%data, problem%exact_lambda)
        if (n < 1) then
            info = -1
            return
        end if
        problem%n = n
        allocate(problem%x(0:n), problem%weights(0:n), &
                 problem%matrix(0:n,0:n), problem%data(0:n), &
                 problem%exact_lambda(0:n), stat=stat)
        if (stat /= 0) then
            info = stat
            call destroy_problem(problem)
            return
        end if

        do i = 0, n
            problem%x(i) = real(i,dp) / real(n,dp)
            problem%weights(i) = 1.0_dp / real(n,dp)
        end do
        problem%weights(0) = 1.0_dp / (2.0_dp*real(n,dp))
        problem%weights(n) = problem%weights(0)

        do j = 0, n
            do i = 0, n
                problem%matrix(i,j) = problem%weights(j) * &
                    exp(problem%x(i)*problem%x(j))
            end do
        end do

        do i = 0, n
            problem%data(i) = (exp(problem%x(i)+1.0_dp)-1.0_dp) / &
                              (problem%x(i)+1.0_dp)
            problem%exact_lambda(i) = exp(problem%x(i))
        end do
    end subroutine build_problem
end module gwr_model

Darauf baut der fachliche Ablauf mit Konditionsanalyse, direkter Lösung und konsistenter Regularisierung auf:

! Derselbe Ablauf wie in MATLAB und Python
type(problem_data) :: problem
real(dp) :: lhs(n+1,n+1), rhs(n+1)
real(dp) :: alpha_direct(n+1), alpha_regularized(n+1)
real(dp) :: condition_a, direct_error, regularized_error

call build_problem(n, problem, info)
call condition_number_2(problem%matrix, condition_a, info)

! Direkte Lösung A alpha = y
call solve_linear_system(problem%matrix, problem%data, &
                         alpha_direct, info)
direct_error = maxval(abs(alpha_direct - problem%exact_lambda))

! Konsistente Regularisierung:
! (beta I + A A) alpha_beta = A y
call build_tikhonov_system(problem%matrix, problem%data, beta, &
                           adjoint_weighted, lhs, rhs, info)
call solve_linear_system(lhs, rhs, alpha_regularized, info)
regularized_error = maxval(abs(alpha_regularized - &
                               problem%exact_lambda))

call destroy_problem(problem)

Die Kommentare im vollständigen Quellcode verweisen auf die jeweils entsprechende Python-Funktion. CMake und CTest prüfen Modellaufbau, stabile Referenzfälle, Regularisierung und beide Adjungierten-Konventionen. Ein unabhängiger Python-Vergleich umfasst 29 Ergebniswerte; die größte relative Differenz beträgt 4.818e-10. Ein zusätzlicher Prüflauf mit Speicher- und Laufzeitdiagnostik blieb ebenfalls fehlerfrei.

Der vollständige Fortran-90/95-Quellcode mit CMake und CTest sowie das gemeinsame Referenz-Repository sind auf GitHub verfügbar.

Als Referenz dienen die Ergebnisse des vorangehenden Sage-Abschnitts. Die neue Gegenüberstellung prüft, ob MATLAB, Python und Fortran diese Referenzfälle reproduzieren.

Numerische Übereinstimmung

Zunächst werden kleine und mittlere Werte von n untersucht. Bei der direkten Lösung sind nur die noch ausreichend stabilen Fälle bis \(n=4\) für einen Vergleich der Implementierungen geeignet. Für feinere Gitter werden die regularisierten Rekonstruktionsfehler gegenüber \(\lambda(t)=\exp(t)\) verglichen.

MATLAB, Python und Fortran reproduzieren die stabilen direkten Fälle und die veröffentlichten regularisierten Sage-Werte. In allen aufgeführten Fällen liefern sie auf die dargestellten fünf signifikanten Stellen dieselben Maximalfehler. Beim feinsten gezeigten Fall mit \(n=128\) und \(\beta=10^{-6}\) liegen auch die ungerundeten Ergebnisse eng beieinander.

Bei der direkten Lösung steigt die Konditionszahl dagegen bereits von ungefähr \(9.1\cdot10^5\) für \(n=4\) auf etwa \(7.9\cdot10^{14}\) für \(n=8\). Ab diesem Bereich können verschiedene BLAS-/LAPACK-Backends deutlich unterschiedliche Lösungsvektoren liefern, obwohl sie dieselbe diskrete Matrix verarbeiten. Diese Abweichung ist ein Symptom der Schlechtgestelltheit und kein belastbarer Geschwindigkeits- oder Qualitätsvergleich zwischen den Programmiersprachen.

Laufzeitvergleich

Für den Laufzeitvergleich werden Matrixaufbau, direkte Lösung, Aufbau des Tikhonov-Systems und regularisierte Lösung getrennt gemessen. Die Gesamtzeit ist der Median der Summe dieser vier Phasen pro Lauf. Nicht enthalten sind Programmstart, Paketinstallation und Datei-Ein-/Ausgabe.

Die veröffentlichten MATLAB-Werte vom 14. August bleiben als historische Referenz stehen. Python wurde am 14. Oktober erneut gemessen und dient als Brücke zwischen den Terminen: Bei \(n=1024\) liegt der frische Python-Wert nur 0.88 % über dem historischen Wert. Fortran ist dort 10.3 % schneller als das frisch gemessene Python und 5.0 % schneller als der historische MATLAB-Wert. Bei kleinen Matrizen dominieren dagegen Aufruf- und Verwaltungsanteile; insbesondere bei \(n=128\) liegt Fortran hinter Python.

Laufzeiten der einzelnen Rechenphasen

Die Phasenmessung zeigt ein differenziertes Bild: Fortran baut die Matrix mit kompilierten Schleifen besonders schnell auf. Bei den drei Operationen der dichten linearen Algebra liegen die Implementierungen enger zusammen, weil ein wesentlicher Teil der Rechenzeit in optimierten BLAS- und LAPACK-Bibliotheken anfällt. Die Messwerte vergleichen deshalb nicht abstrakt Programmiersprachen, sondern stets die getesteten Programme zusammen mit Compiler und numerischem Backend.

Reproduzierbarkeit

• identische Problemdefinition und identische Eingabedaten,
• Matrixdimension \(N=n+1\) ausdrücklich angegeben,
• MATLAB-, Python-, NumPy-, Fortran-Compiler- und Sage-Version dokumentiert,
• BLAS-/LAPACK-Backend dokumentiert,
• Threadzahl und CPU-Affinität kontrolliert,
• Warm-up vor der Messung,
• mehrere Wiederholungen mit Median und Streuung,
• Validierungs- und Performancefälle getrennt ausgewiesen.

MATLAB wurde mit -singleCompThread gestartet. Für Python und Fortran waren OPENBLAS_NUM_THREADS=1, OMP_NUM_THREADS=1 und MKL_NUM_THREADS=1 gesetzt. Alle Prozesse liefen mit taskset -c 0 auf demselben logischen CPU-Kern.

Der Einthreadmodus wurde bewusst als klar definierte Vergleichsbasis gewählt. Dadurch dominieren unterschiedliche automatische Threadstrategien von MKL und OpenBLAS die Messung nicht. Eine separate Messreihe untersucht im Folgenden die Parallelskalierung mit identischer Threadzahl und kontrollierter CPU-Affinität.

Threadskalierung auf physischen Kernen

Für die Threadmessung wird die dichte Performanceaufgabe auf \(n=8192\), also \(N=8193\), vergrößert. Gemessen werden 1, 2, 4, 8, 16, 20 und 40 Rechenthreads. Für MATLAB und Python liegen jeweils fünf unabhängige Läufe vom 14. August vor; Fortran wurde am 14. Oktober in drei unabhängigen Läufen gemessen. Bis 20 Threads liegen alle verwendeten physischen Kerne auf dem ersten CPU-Sockel. Der 40-Thread-Punkt verwendet die 40 physischen Kerne beider Sockel. Hyperthreading bleibt vollständig ausgeschlossen.

Die Parallelisierung wird nicht im Modellcode mit parfor, multiprocessing oder OpenMP programmiert. MATLAB begrenzt seine Rechenthreads mit maxNumCompThreads; bei Python und Fortran steuert OPENBLAS_NUM_THREADS das OpenBLAS-Backend. taskset bindet die Prozesse an dieselben physischen CPUs. Matrixprodukte und lineare Löser werden damit innerhalb von MKL beziehungsweise OpenBLAS parallel ausgeführt.

MATLAB, Beispiel mit 16 Threads:
taskset -c 0-15 matlab -batch "maxNumCompThreads(16); ..."

Python/NumPy, Beispiel mit 16 Threads:
OPENBLAS_NUM_THREADS=16 OMP_NUM_THREADS=16 MKL_NUM_THREADS=16 \
taskset -c 0-15 python -m gwr_parameteridentifikation.benchmark ...

Fortran/OpenBLAS, Beispiel mit 16 Threads:
OPENBLAS_NUM_THREADS=16 OMP_NUM_THREADS=16 \
taskset -c 0-15 ./gwr_thread_benchmark_16

Ergebnisse der Threadskalierung

Alle drei Implementierungen erreichen bis acht Threads ungefähr den fünffachen Durchsatz. Fortran erzielt mit 16 Threads eine Gesamtzeit von 3.992 s und damit einen Speedup von 5.64. Dies ist zugleich der beste Fortran-Messpunkt; bei 20 Threads steigt die Zeit geringfügig auf 4.130 s.

Die Fortran-Phasenmessung zeigt die Wirkung des aktuellen Backends: Bei 16 Threads beschleunigt sich der mit DGEMM berechnete Aufbau des Tikhonov-Systems um den Faktor 7.71. Die direkte und die regularisierte DGESV-Lösung erreichen Speedups von 4.53 beziehungsweise 4.76. Der explizite Matrixaufbau bleibt einthreadig, macht mit rund 0.23 s jedoch nur einen kleinen Teil der Gesamtlaufzeit aus.

Eine Kontrollmessung mit dem zuvor installierten OpenBLAS 0.3.29 hatte für DGESV nahezu konstante Zeiten ergeben. Ursache war ein inzwischen korrigierter Schwellwertfehler dieser Version: Bei einer rechten Seite wurde die Parallelisierung anhand von N × NRHS statt anhand der quadratischen Matrixgröße aktiviert. Mit OpenBLAS 0.3.34 skaliert derselbe unveränderte Fortran-95-Code erwartungsgemäß.

Auch der Matrixaufbau reagiert unterschiedlich: MATLAB parallelisiert den vektorisierten Aufbau in dieser Konfiguration sichtbar, während NumPys outer/exp-Folge und die expliziten Fortran-Schleifen einthreadig bleiben. Der 40-Kern-Fall lief unter der Standard-NUMA-Policy von Linux und ist als eigener Zwei-Sockel-Messpunkt zu lesen. Die Skalierungskurve bestätigt damit, dass die konkrete Backend- und Speicherarchitektur wichtiger ist als die bloße Angabe einer maximalen Threadzahl.

Was der Vergleich zeigt

MATLAB, Python und Fortran können dieselbe Parameteridentifikation mit praktisch übereinstimmenden Resultaten abbilden. Entscheidend für die fachliche Qualität sind jedoch eine konsistente Diskretisierung, eine geeignete Regularisierung und eine sachgerechte Beurteilung der Kondition des diskretisierten Problems. Die Wahl der Programmiersprache beeinflusst Bedienung und konkrete Laufzeit, ändert aber nichts an der mathematischen Schlechtgestelltheit.

Bei großen Diskretisierungen zeigt sich außerdem, dass die Wahl des Algorithmus wesentlich mehr bewirken kann als der Wechsel zwischen MATLAB, Python und Fortran. Der Kern \(\exp(x_i x_j)\) besitzt eine separierbare Entwicklung. Dadurch lässt sich die vollständige Kernmatrix durch voll besetzte Faktormatrizen mit nur wenigen Spalten darstellen.

Ausblick: Niedrigrangrechnung ohne vollständige Matrix

Die daraus entstehende Niedrigrangrechnung mit separierbaren Kernen ist ein eigenständiges Thema. Im vorliegenden Beispiel ersetzen voll besetzte Faktormatrizen mit nur wenigen Spalten die vollständige Kernmatrix. Damit folgt die Rechnung derselben strukturorientierten Grundidee, die auch hierarchischen Matrizen und dem Panel Clustering zugrunde liegt: Die Darstellung speichert nicht alle Matrixeinträge einzeln, sondern bildet die Operatorwirkung aus kompakten Faktoren.

Hinweis: Die ergänzten Abschnitte zu Sage, MATLAB, Python und Fortran sowie zur Threadskalierung wurden mit Unterstützung KI-gestützter Werkzeuge ausgearbeitet. Mathematische Aussagen, Implementierungen und Messergebnisse wurden fachlich geprüft und durch reproduzierbare Rechnungen validiert.