std.decomp — Matrixzerlegungen

Echte SVD nach einseitigem Jacobi, Pseudoinverse, Rang-k-Näherung, polare Zerlegung und allgemeine Eigenwerte über QR-Iteration.

→ Mathematik-Units · std.linalg · std.sparse

import std.decomp;

Schicht 1, hängt an std.linalg, std.math und std.alloc. 14 Funktionen — die kleinste Schnittstelle der Sammlung bei 913 Zeilen Umsetzung. Alle Beispiele mit lyxc 1.0.21A übersetzt und ausgeführt.

Verhältnis zu ''std.linalg''

std.decomp ist die Fortsetzung von std.linalg, nicht dessen Ersatz. Die Trennung verläuft entlang „Lösen„ gegen „Zerlegen“:

Aufgabe	Unit
Quadratisches System lösen, Determinante, Inverse	std.linalg
Symmetrische Eigenwerte (Jacobi)	`std.linalg`
SVD, Pseudoinverse, Rang-k-Näherung	`std.decomp`
Allgemeine (nicht symmetrische) Eigenwerte	`std.decomp`

Beides zusammen wären über 2000 Zeilen in einer Unit. Getrennt bleibt die Aufteilung verständlich, und wer nur Gleichungssysteme löst, zieht die schwerere Unit nicht mit.

SVD und was daraus folgt

unit main;
import std.io;
import std.string;
import std.alloc;
import std.linalg;
import std.decomp;

fn ZF(t: pchar, v: f64): void { Print(t); PrintF64(v); }
fn ZI(t: pchar, v: int64): void { Print(t); PrintLn(IntToStr(v)); }

fn main(): int64 {
    var m: int64 := 3; var n: int64 := 2;
    var a: int64 := alloc(m * n * 8);
    LinSet(a, n, 0,0, 1.0); LinSet(a, n, 0,1, 2.0);
    LinSet(a, n, 1,0, 3.0); LinSet(a, n, 1,1, 4.0);
    LinSet(a, n, 2,0, 5.0); LinSet(a, n, 2,1, 6.0);

    var s: int64 := alloc(n * 8);
    DecompSvdValues(a, m, n, s);
    ZF("Singulaerwert 0 : ", peekf64(s));
    ZF("Singulaerwert 1 : ", peekf64(s + 8));
    ZF("Norm2           : ", DecompNorm2(a, m, n));
    ZF("Cond2           : ", DecompCond2(a, m, n));
    ZI("Rank            : ", DecompRank(a, m, n, DecompEpsilon()));

    var b: int64 := alloc(m * 8);
    pokef64(b, 1.0); pokef64(b + 8, 2.0); pokef64(b + 16, 3.0);
    var x: int64 := alloc(n * 8);
    var ok: bool := DecompLstSq(a, m, n, b, DecompEpsilon(), x);
    Print("LstSq ok="); Print(IntToStr(ok as int64)); PrintLn("");
    ZF("x0              : ", peekf64(x));
    ZF("x1              : ", peekf64(x + 8));

    var c: int64 := alloc(4 * 8);
    LinSet(c, 2, 0,0, 4.0); LinSet(c, 2, 0,1, 1.0);
    LinSet(c, 2, 1,0, 2.0); LinSet(c, 2, 1,1, 3.0);
    ZF("Spektralradius  : ", DecompSpectralRadius(c, 2));

    var evr: int64 := alloc(2 * 8);
    var evi: int64 := alloc(2 * 8);
    DecompEigenGeneral(c, 2, evr, evi);
    ZF("Eigenwert 0 Re  : ", peekf64(evr));
    ZF("Eigenwert 1 Re  : ", peekf64(evr + 8));
    return 0;
}

Singulaerwert 0 : 9.525518
Singulaerwert 1 : 0.514300
Norm2           : 9.525518
Cond2           : 18.521305
Rank            : 2
LstSq ok=1
x0              : 0.000000
x1              : 0.499999
Spektralradius  : 5.000000
Eigenwert 0 Re  : 5.000000
Eigenwert 1 Re  : 2.000000

Sämtliche Werte stimmen mit der Referenzrechnung überein.

Drei Größen fallen mit den Singulärwerten ab, ohne dass man sie einzeln rechnen müsste:

DecompNorm2 ist der größte Singulärwert (9,525518).
DecompCond2 ist das Verhältnis von größtem zu kleinstem (9,525518 / 0,514300 = 18,52).
DecompRank ist die Anzahl der Singulärwerte oberhalb der Toleranz.

Rang ist eine Frage der Toleranz

<WRAP alert> DecompRank braucht ein tol, und die Wahl ist eine Entscheidung. In exakter Arithmetik ist der Rang eindeutig; in f64 ist ein Singulärwert von 1e−16 nicht von null zu unterscheiden. DecompEpsilon() liefert 1e−15 als Vorgabe.

Bei gemessenen Daten ist die richtige Schwelle größer: dort trennt sie Signal von Rauschen, und ein Singulärwert unterhalb der Messgenauigkeit ist keine echte Dimension. Wer die Vorgabe nimmt, bekommt bei verrauschten Daten regelmäßig vollen Rang — mathematisch korrekt, sachlich nutzlos. </WRAP>

Pseudoinverse und Ausgleichsrechnung

Funktion	Zweck
`DecompPseudoInverse(a, m, n, tol, out)`	Moore-Penrose-Inverse
`DecompLstSq(a, m, n, b, tol, x)`	überbestimmtes System, kleinste Quadrate
`DecompLowRank(a, m, n, k, out)`	beste Näherung vom Rang `k`

Die Pseudoinverse gibt es auch für nicht quadratische und singuläre Matrizen, wo LinInverse scheitert. Für ein überbestimmtes System liefert DecompLstSq die Lösung mit dem kleinsten Fehlerquadrat — oben x = (0, 0.5) für drei Gleichungen mit zwei Unbekannten.

out darf nicht auf a zeigen. DecompPseudoInverse und DecompPolar prüfen das und geben false zurück; verlassen sollte man sich darauf nicht, sondern getrennte Puffer verwenden. DecompPseudoInverse setzt außerdem m >= n voraus.

DecompLowRank ist die Grundlage von Kompression und Entrauschen über die SVD: die k größten Singulärwerte tragen den Hauptteil der Information, der Rest ist oft Rauschen.

Eigenwerte allgemeiner Matrizen

Funktion	Zweck
`DecompHessenberg(a, n, h, q)`	Vorbereitung für die QR-Iteration
`DecompEigenGeneral(a, n, evRe, evIm)`	Eigenwerte, komplex
`DecompEigenVector(a, n, lambda, vec)`	Eigenvektor zu einem bekannten Eigenwert
`DecompSpectralRadius(a, n)`	größter Betrag aller Eigenwerte
`DecompPolar(a, n, q, p)`	polare Zerlegung `A = Q·P`

Die Eigenwerte kommen in zwei Feldern — Real- und Imaginärteil. Eine nicht symmetrische reelle Matrix kann komplexe Eigenwerte haben, und zwar immer paarweise konjugiert. Wer nur evRe liest, bekommt bei einer Drehmatrix zwei gleiche Zahlen und übersieht, dass die eigentliche Information im Imaginärteil steckt.

Für symmetrische Matrizen ist LinEigenSym aus std.linalg die bessere Wahl: dort sind alle Eigenwerte reell, das Jacobi-Verfahren ist stabiler und liefert die Eigenvektoren gleich mit.

DecompSpectralRadius beantwortet die Frage, ob eine Iteration x ← A·x konvergiert: unter 1 ja, darüber nicht.

Letzte Aktualisierung: 2026-08-16 · alle Beispiele mit lyxc 1.0.21A übersetzt und ausgeführt; Singulärwerte, Kondition, Rang, Ausgleichslösung und Eigenwerte gegen eine unabhängige Referenzrechnung geprüft.