[java] Moindres carrés pour systèmes lineaires

Version imprimable

[java] Fitting de courbe (moindres carrés)

Le "fitting" de courbe (ou "ajustement" en francais) a pour objectif de trouver la meilleure courbe passant par un ensemble de points.

Par "meilleure courbe" on entend celle qui minimise une erreur. Dans notre cas, cette erreur est l'ecart quadratique entre la courbe et les points.

Lorsque la courbe recherchée est une forme linéaire, la methode des "moindres carrés" permet de calculer les coefficients de la courbe en résolvant un système linéraire.

Le code:

La classe Matrix, qui s'occupe de stocker et inverser une matrice:

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
 
/**
 * Structure et Operations sur une matrice
 * 
 * @author Xavier PHILIPPEAU
 *  
 */
public class Matrix {
 
	// les données de la matrice ordonnnée en [x=colonne][y=ligne]
	public double value[][];
 
	// les dimensions de la matrice
	public final int rows, cols;
 
	// constructeur de matrice quelconque (vide)
	public Matrix(int cols, int rows) {
		this.cols = cols;
		this.rows = rows;
		this.value = new double[this.cols][this.rows];
	}
 
	// constructeur de matrice carré (pré-remplies)
	public Matrix(int size, double[][] squarearray) {
		this.cols = size;
		this.rows = size;
		this.value = squarearray;
	}
 
	/**
         * inversion par pivot de Gauss
         *
         * @return l'inverse de la matrice (ou null si non inversible)
         */
	public Matrix inverse() {
 
		//Une matrice ne peut etre inversée que si elle est carrée.
		if (this.rows != this.cols) {
			System.err.println("Inversion : La matrice n'est pas carré !!");
			return null;
		}
 
		// Création de la matrice de travail T = [ this | Identité ]
		Matrix T = new Matrix(this.cols * 2, this.cols);
		for (int y = 0; y < this.rows; y++) {
			for (int x = 0; x < this.cols; x++) {
				T.value[x][y] = this.value[x][y];
				if (x == y)	T.value[this.cols + x][y] = 1;
			}
		}
 
		// Pour chaque ligne de la matrice T
		for (int x = 0; x < T.rows; x++) {
 
			// Cherche la ligne avec le pivot max (en valeur absolue)
			int bestline = x;
			double pivot = T.value[x][bestline];
			for (int y = x + 1; y < T.rows; y++) {
				if (Math.abs(T.value[x][y]) > Math.abs(pivot)) {
					bestline = y;
					pivot = T.value[x][bestline];
				}
			}
 
			if (pivot == 0) {
				System.err.println("Inversion : Le pivot est nul,inversion impossible !!");
				return null;
			}
 
			// Echange des lignes (si necessaire)
			if (bestline != x) {
				double tmp;
				for (int t = 0; t < T.cols; t++) {
					tmp = T.value[t][x]; T.value[t][x] = T.value[t][bestline]; T.value[t][bestline] = tmp;
				}
			}
 
			// Normalisation de la ligne du pivot
			for (int t = 0; t < T.cols; t++) T.value[t][x] /= pivot;
 
			// elimination des autres lignes
			for (int y = 0; y < T.rows; y++) {
				if (y==x) continue;
				double coef = T.value[x][y];
				for (int t = 0; t < T.cols; t++) T.value[t][y] -= coef * T.value[t][x];
			}
		}
 
		// recupere la partie droite de T qui contient l'inverse de la matrice
		// (la partie gauche devrait contenir l'identité)
		Matrix inverse = new Matrix(this.cols,this.rows);
		for (int y = 0; y < this.rows; y++)
			for (int x = 0; x < this.cols; x++)
				inverse.value[x][y] = T.value[this.cols+x][y];
 
		return inverse;
	}
 
	/**
         * Multiplication de la matrice par un vecteur
         *
         * @param X le vecteur
         * @return le vecteur résultat: this[][] * X[]
         */
	public double[] multiply(double[] X) {
		double[] Y = new double[this.rows];
		for (int y = 0; y < this.rows; y++)
			for (int x = 0; x < this.cols; x++)
				Y[y] += this.value[x][y]*X[x];
		return Y;
	}	
 
}

La classe LinearEquation qui s'occupe de modeliser une equation lineaire (coefs + valeur)

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
 
/**
 * Structure d'une equation lineaire
 * 
 *  coef[0]*X0 +  coef[1]*X1 + ... + coef[n]*Xn = residu
 *  
 *  ou X=(X0,X1,...,Xn) est le vecteur "inconnu"
 *
 * @author Xavier PHILIPPEAU
 */
public class LinearEquation {
 
	/**
         * Les coefficients de l'equation
         */
	public double[] coefficents;
 
 
	/**
         * Le résidu (= le résultat) de l'equation 
         */
	public double residual;
 
	/**
         * Constructeur
         * 
         * @param unknownsCount Le nombre d'inconnues
         */
	public LinearEquation(int unknownsCount) {
		this.coefficents = new double[unknownsCount];
	}
 
	/**
         * Constructeur
         * 
         * @param coefficents Les coefficients de l'equation
         * @param residual Le résidu (= le résultat) de l'equation 
         */
	public LinearEquation(double[] coefficents, double residual) {
		this.coefficents = coefficents;
		this.residual = residual;
	}
 
}

La classe LeastSquare qui s'occupe de construire et résoudre le système d'equations lineaires par la méthode des "moindres carrés":

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
 
import java.util.List;
 
/**
 * résolution du système d'equations par la methode des moindres carrés
 *
 * @author Xavier PHILIPPEAU
 */
public class LeastSquare {
 
	/**
         * Resolution d'un systeme lineaire
         * 
         * @param systemsize Nombre d'inconnues dans le systeme
         * @param system Le systeme d'equations lineaires
         * @return le vecteur solution du système
         */
	public static double[] linearFit(int systemsize, List<LinearEquation> system) {
 
		// taille du systeme a resoudre (= nombre d'inconnues)
		int n = systemsize;
 
		// On cherche a trouver X qui resoud le systeme:
		//
		// J.X = R
		//
		// on multiplie par Jt, la transposée de J:
		//
		// Jt.J.X = Jt.R
		//
		// on pose:
		//                    n
		// A[c][l] = Jt.J = Somme[ coefficient[c][k] * coefficient[l][k] ]
		//                   k=0
		//
		//                 n
		// Y[l] = Jt.R = Somme[ residual[k] * coefficient[l][k] ]
		//                k=0
		//
		// et "n", le nombre d'equation dans le systeme.
		//
		// Le systeme a resoudre devient ainsi:
		//
		// A.X = Y
 
		Matrix A = new Matrix(n,n);
		double Y[] = new double[n];
		for(LinearEquation equation:system) {
			for (int l = 0; l < n; l++) {
				for (int c = 0; c < n; c++) {
					A.value[c][l] += equation.coefficents[c] * equation.coefficents[l];
				}
				Y[l] += equation.residual * equation.coefficents[l];
			}
		}
 
		// Resolution par inversion:
		//
		// A . X = Y  ==>  X = A^-1 . Y
 
		double X[] = (A.inverse()).multiply(Y);		
 
		return X;
	}
 
}

PRomu@ld : sources

Exemples d'utilisation

Trouver la "meilleure ligne" au sens des moindres carrés:

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
 
public static double[] bestLine(List<double[]> data) {
 
	List<LinearEquation> system = new ArrayList<LinearEquation>();
 
	for(double[] d : data) {
		double x = d[0];
		double y = d[1];
 
		// a.1 + b.X = Y
		//
		// (X,Y) sont connus (données du probleme)
		// (a,b) sont les 2 inconnues
 
		LinearEquation equation = new LinearEquation(2);
		equation.coefficents = new double[] { 1, x };
		equation.residual = y;
 
		system.add(equation);
	}
 
	// resolution du systeme lineaire
	double[] C = LeastSquare.linearFit(2, system);
	return C;
 
	// La solution s'ecrit: Y = C[0] + C[1]*X
}

Trouver la "meilleure parabole" au sens des moindres carrés:

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
 
public static double[] BestParabol(List<double[]> data) {
 
	List<LinearEquation> system = new ArrayList<LinearEquation>();
 
	for(double[] d : data) {
		double x = d[0];
		double y = d[1];
 
		// a.1 + b.X + c.X^2 = Y
		//
		// (X,Y) sont connus (données du probleme)
		// (a,b,c) sont les 3 inconnues
 
		LinearEquation equation = new LinearEquation(3);
		equation.coefficents = new double[] { 1, x, x*x };
		equation.residual = y;
 
		system.add(equation);
	}
 
	// resolution du systeme lineaire
	double[] C = LeastSquare.linearFit(3, system);
	return C;
 
	// La solution s'ecrit: Y = C[0] + C[1]*X + C[2]*X^2
}

Trouver le "meilleure polynome de degré n" au sens des moindres carrés:

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
 
public static double[] BestPolynom(int order, List<double[]> data) {
 
	List<LinearEquation> system = new ArrayList<LinearEquation>();
 
	for(double[] d : data) {
		double x = d[0];
		double y = d[1];
 
		// a.1 + b.X + c.X^2 + ... + z.X^ordre = Y
		//
		// (X,Y) sont connus (données du probleme)
		// (a,b,c,...z) sont les inconnues
 
		LinearEquation equation = new LinearEquation(order+1);
 
		equation.residual = y;
		for(int n=0;n<=order;n++)
			equation.coefficents[n] = Math.pow(x, n);
 
		system.add(equation);
	}
 
	// resolution du systeme lineaire
	double[] C = LeastSquare.linearFit(order+1, system);
	return C;
 
	// La solution s'ecrit: Y = C[0] + C[1]*X + C[2]*X^2 + ... + C[n]*X^n
}

un exemple/explication du cas "non lineaire" appliqué à un cercle:
http://www.orbitals.com/self/least/least.htm

10/07/2007, 12h10
souviron34

;) je la traduirai en C (Java moi yen a pas connaître) et je posterais le texte du papier cette semaine (peut-être aujourdhui ?) .. (je viens de rentrer de voyage..)

Trouver la "meilleure ellipse" au sens des moindres carrés, en utilisant les formules de mathworld.wolfram.com:

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
 
public static double[] BestEllipsoid(List<double[]> data) {
 
	List<LinearEquation> system = new ArrayList<LinearEquation>();
 
	for(double[] d : data) {
		double x = d[0];
		double y = d[1];
 
		// a.x^2 + 2b.xy + c.y^2 + 2d.x + 2f.y + g = 0
		//
		// (X,Y) sont connus (données du probleme)
		// on pose g = -1
		// (a,b,c,d,f) sont les 5 inconnues
 
		LinearEquation equation = new LinearEquation(5);
		equation.coefficents = new double[] { x*x, x*y, y*y, x, y };
		equation.residual = 1.0;
 
		system.add(equation);
	}
 
	// resolution du systeme lineaire
	double[] C = LeastSquare.linearFit(5, system);
 
	// evaluation des 5 inconnues
	double a = C[0];
	double b = C[1]/2;
	double c = C[2];
	double d = C[3]/2;
	double f = C[4]/2;
	double g = -1;
 
	// calculs des caracteristiques de l'ellipse:
	//
	// d'après http://mathworld.wolfram.com/Ellipse.html
	//
	// formules (15), (16), ..., (23)
 
	// verifie que l'ellipse est valide
	double delta=a*(c*g-f*f)-b*(b*g-d*f)+d*(b*f-d*c);
	double J=a*c-b*b;
	double I=a+c;
	if(delta==0) return null;
	if(J<=0) return null;
	if((delta/I)>=0) return null;
 
	/* calcul du centre de l'ellipse */
	double xcentre = (c*d-b*f)/(b*b-a*c);
	double ycentre = (a*f-b*d)/(b*b-a*c);
 
	/* angle de rotation (en radian) */
	double angle = 0.5*Math.atan(2*b/(c-a));
 
	/* longueurs des 2 demi-axes */
	double num = 2*(a*f*f + c*d*d + g*b*b - 2*b*d*f - a*c*g);
	double sq = Math.sqrt( 1 + 4*b*b/((a-c)*(a-c)) );
 
	double denom1 = (b*b-a*c) * ( (c-a) * sq - (c+a) );
	double xhalflength = Math.sqrt(num/denom1);
 
	double denom2 = (b*b-a*c) * ( (a-c) * sq - (c+a) );
	double yhalflength = Math.sqrt(num/denom2);
 
	/* retourne toujours un angle positif  */
	if (angle<0) {
		angle+=Math.PI/2;
		// echange les longueurs des demi-axes
		double tmp=xhalflength;
		xhalflength=yhalflength;
		yhalflength=tmp;
	}
 
	return new double[] { xcentre, ycentre, angle, xhalflength, yhalflength};
}

Et un petit test pour s'assurer que ca marche...

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
 
// construction des points d'une ellipse
List<double[]> data = new ArrayList<double[]>();
for(int i=0;i<100;i++) {
	double t = i/100.0;
	double a = t*2*Math.PI;
 
	// ellipse de demi-axe 17 x 10, centrée en (0,0)
	double xe = 17*Math.cos(a);
	double ye = 10*Math.sin(a);
 
	// + une rotation de 45 degrés
	double r = Math.toRadians(45);
	double xr = xe*Math.cos(r)+ye*Math.sin(r);
	double yr = -xe*Math.sin(r)+ye*Math.cos(r);
 
	// + une translation de (15,8)
	double x = xr+15;
	double y = yr+8;
 
	data.add( new double[] {x,y});
}
 
double[] E = BestEllipsoid(data);
 
System.out.println("centre: "+E[0]+","+E[1]);
System.out.println("angle: "+Math.toDegrees(E[2])+" (degrés)");
System.out.println("demi-axe: "+E[3]+" x "+E[4]);
 
// affiche:
// centre: 14.999999999999954,8.000000000000052
// angle: 44.99999999999916 (degrés)
// demi-axe: 17.00000000000003 x 10.000000000000034

17/09/2007, 19h27
FR119492

[java] Fitting de courbe (moindres carrés)

Deux remarques à l'intention de pseudocode:

1) Si la tâche à accomplir consiste à résoudre un système linéaire, ce qui est très souvent le cas, il est inutilement coûteux en temps d'exécution d'inverser la matrice (bien que ce soit ce que font la plupart des gens); le plus efficace est de la factoriser par la méthode LU dans le cas général, ou par la méthode de Cholesky, deux fois plus rapide, si elle est symétrique définie positive.

2) Si tu cherches le vecteur qui satisfait le moins mal possible (au sens des moindres carrés) un système linéaire surdéterminé, regarde la méthode des valeurs singulières (singular value decomposition); tu trouveras toute la documentation à ce sujet dans Numerical Recipes ou dans la description de la fonction svd de MatLab.

A ta disposition si tu as encore des questions à ce sujet.
13/02/2008, 00h28
souheyeb

bonjour,
je souhaite demander si je peux utiliser cette méthode pour résoudre un système de la forme : X^t*A=o , avec A matrice non carrée(n*m)
merci d'avance
13/02/2008, 09h44
pseudocode

Citation:

Envoyé par souheyeb

bonjour,
je veus demander si je peut utiliser cette methode pour resoudre un systeme de la forme : X^t*A=o , avec A matrice non carrée(n*m)
merci d'avance

Heu... Non. La méthode des moindres carrés a pour objectif de trouver les "meilleurs" coefficients d'une courbe. La résolution du système linéaire surdéterminé est l'une des étapes de cette méthode, mais nullement le but en soi. D'ailleurs, comme l'a fait remarquer FR119492, elle n'est meme pas obligatoire.

Pour la résolution d'un système linéaire surdéterminé, je te conseille plutôt la factorisation QR ou la décomposition SVD.
10/05/2012, 16h55
pservolex

Pour passer en 3d et chercher l'equation d'un plan , il suffit de déclarer
double z = d[2];

et
de chercher l'equation lineaire z=aX+bY+c ?

UN petit coup de main pour un informaticien non mathematicien ? (pitié!)

Une librairie ?

Bonjour,

Merci pour ce sujet, pseudocode.
On peut aussi utiliser une librairie d'algèbre linéaire pour créer le système linéaire et le résoudre... (ça fait moins de code à écrire et il y a peut-être des améliorations pour la résolution du système quand il est mal conditionné, par rapport à tA.A x = tA.b )
Jama, par exemple ?
http://math.nist.gov/javanumerics/jama/

Exemple :
Code:

1 2 3 4 5 6 7 double[][] array = {{1.,2.,3},{4.,5.,6.},{7.,8.,10.}}; Matrix A = new Matrix(array); Matrix b = Matrix.random(3,1); Matrix x = A.solve(b); Matrix Residual = A.times(x).minus(b); double rnorm = Residual.normInf();

Version avec Jama

Citation:

Envoyé par FR119492

1) Si la tâche à accomplir consiste à résoudre un système linéaire, ce qui est très souvent le cas, il est inutilement coûteux en temps d'exécution d'inverser la matrice (bien que ce soit ce que font la plupart des gens); le plus efficace est de la factoriser par la méthode LU dans le cas général, ou par la méthode de Cholesky, deux fois plus rapide, si elle est symétrique définie positive.

A priori, pour une résolution "au mieux" (matrice rectangulaire, nb lignes > nb colonnes), Jama utilise une décomposition QR...

Code:

1
2
3
4
5
 
   public Matrix solve (Matrix B) {
      return (m == n ? (new LUDecomposition(this)).solve(B) :
                       (new QRDecomposition(this)).solve(B));
   }

J'ai écrit une version dérivée du code de Xavier, adaptée à l'utilisation de Jama, voilà ce que ça donne :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
 
import Jama.Matrix;
import java.awt.geom.Point2D;
import java.util.AbstractList;
 
[...]
 
	/**
         * Renvoie les paramètre de l'ellipse qui s'ajuste mieux au nuage de
         *  points 2D reçu (au sens des moindres carrés)
         * @param lstPts2D liste de points 2D (Point2D.Double) du contour de l'ellipse
         * @return un tableau qui contient les paramètres de l'ellipse
         */
	public static double[] getParametres_ellipse_mc( AbstractList lstPts2D ) {
		// construction du système linéaire
		int nNb_pts = lstPts2D.size();
		Matrix mA = new Matrix( nNb_pts, 5 );
		Matrix vB = new Matrix( nNb_pts, 1 );
		for( int i = 0; i < nNb_pts; ++i ) {
			Point2D.Double pt = (Point2D.Double)lstPts2D.get( i );
			mA.set( i, 0, pt.x * pt.x );
			mA.set( i, 1, pt.x * pt.y );
			mA.set( i, 2, pt.y * pt.y );
			mA.set( i, 3, pt.x );
			mA.set( i, 4, pt.y );
			vB.set(i, 0, -1.0 );
		}
		// resolution du systeme lineaire A.x = b
		Matrix vX = mA.solve( vB );
 
		// evaluation des 5 inconnues
		double a = vX.get( 0, 0 );
		double b = vX.get( 1, 0 ) / 2.0;
		double c = vX.get( 2, 0 );
		double d = vX.get( 3, 0 ) / 2.0;
		double f = vX.get( 4, 0 ) / 2.0;
		double g = 1.0;
 
		// calcule les caracteristiques de l'ellipse:
		//
		// d'après http://mathworld.wolfram.com/Ellipse.html
		//
		// formules (15), (16), ..., (23)
 
		// verifie que l'ellipse est valide
		double delta = a * ( c * g - f * f) - b * ( b * g - d * f ) + d * ( b * f - d * c );
		double nJ = a * c - b * b;
		double nI = a + c;
		if ( delta == 0 ) return null;
		if ( nJ <= 0 ) return null;
		if ( ( delta / nI ) >= 0) return null;
 
		// calcul du centre de l'ellipse
		double nX_centre = ( c*d - b*f )/( b*b - a*c );
		double nY_centre = ( a*f - b*d )/( b*b - a*c );
 
		// angle de rotation (en radian)
		double nAngle_rad = 0.5 * Math.atan( 2.0 * b / ( c-a ));
 
		// longueurs des 2 demi-axes
		double num = 2.0 * ( a * f*f + c * d*d + g* b*b - 2*b*d*f - a*c*g);
		double sq = Math.sqrt( (a-c)*(a-c) + 4.0*b*b );
 
		double denom1 = ( b*b - a*c ) * ( sq - ( a + c ) );
		double nDemi_axe_x = Math.sqrt( num / denom1 );
 
		double denom2 = ( b*b - a*c ) * ( -1.0 * sq - ( a + c ) );
		double nDemi_axe_y = Math.sqrt( num / denom2 );
 
		// retourne un angle entre [0 et 2\Pi[
		if ( nAngle_rad < 0.0 ) {
			nAngle_rad += 2.0 * Math.PI;
		}
		return new double[]{nX_centre, nY_centre, nAngle_rad, nDemi_axe_x, nDemi_axe_y };
	}// getParametres_ellipse_mc

Avec la même procédure de test que Xavier j'obtiens aussi les bons résultats :

Code:

1
2
3
4
 
centre: 14.999999999999991,7.999999999999996
angle: 44.99999999999999 (degrés)
demi-axe: 17.000000000000025 x 9.999999999999982

Merci encore pour ce sujet, il m'a bien aidé...

20/08/2012, 11h14
Invité

La version "EllipseFitting" de Michael Doube est intéressante (voir le code ici) :
http://github.com/mdoube/BoneJ/blob/...itEllipse.java

Voir notamment la fonction "varToDimensions"...
20/08/2012, 11h33
souviron34

tu peux regarder j'en ai mise une ici-même en Fortran, facilement convertible.. :)
20/08/2012, 11h59
Aleph69

Bonjour,

inverser une matrice pour ensuite faire

Citation:

// Resolution par inversion:
//
// A . X = Y ==> X = A^-1 . Y

double X[] = (A.inverse()).multiply(Y);

c'est vraiment très maladroit!
23/09/2012, 14h39
saturn1

C'est exactement ce que je cherchais.

(enfin je cherchais un plane fitting)

Merci.
24/04/2013, 22h31
ToTo13

En revanche, je pense que cette méthode doit pouvoir être utilisée pour trouver la meilleure gaussienne 2D passant par un ensemble de points.
Ou ai je tort ?

[EDIT] : J'ai définitivement tort car l'équation de la Gaussienne n'est pas linéaire à cause de l'exponentielle.
Donc comment modifier la méthode pour résoudre dans le cas Gaussien ?