UpsilonNumworks/Upsilon
1
0
Fork 0
mirror of https://github.com/UpsilonNumworks/Upsilon.git synced 2026-03-19 05:40:38 +01:00
Code Issues Packages Projects Releases Wiki Activity

[apps/reg] Cap the number of matrix modifications

Browse Source
This commit is contained in:
Léa Saviot
2018-06-18 17:48:59 +02:00
committed by Émilie Feral
parent c56652cda9
commit f031a642b8
2 changed files with 16 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files

19
apps/regression/model/model.cpp
View File

@@ -65,7 +65,9 @@ void Model::fitLevenbergMarquardt(Store * store, int series, double * modelCoeff
operandsB[j] = betaCoefficient(store, series, modelCoefficients, j);
}
double modelCoefficientSteps[Model::k_maxNumberOfCoefficients];
solveLinearSystem(modelCoefficientSteps, coefficientsAPrime, operandsB, n, context);
if (solveLinearSystem(modelCoefficientSteps, coefficientsAPrime, operandsB, n, context) < 0) {
break;
}
// Compare new chi2 with the previous value
double newModelCoefficients[Model::k_maxNumberOfCoefficients];
for (int i = 0; i < n; i++) {
@@ -146,7 +148,7 @@ double Model::betaCoefficient(Store * store, int series, double * modelCoefficie
return result;
}
void Model::solveLinearSystem(double * solutions, double * coefficients, double * constants, int solutionDimension, Context * context) {
int Model::solveLinearSystem(double * solutions, double * coefficients, double * constants, int solutionDimension, Context * context) {
int n = solutionDimension;
assert(n <= k_maxNumberOfCoefficients);
double coefficientsSave[k_maxNumberOfCoefficients * k_maxNumberOfCoefficients];
@@ -154,18 +156,25 @@ void Model::solveLinearSystem(double * solutions, double * coefficients, double
coefficientsSave[i] = coefficients[i];
}
int inverseResult = Matrix::ArrayInverse(coefficients, n, n);
if (inverseResult < 0) {
int numberOfMatrixModifications = 0;
while (inverseResult < 0 && numberOfMatrixModifications < k_maxMatrixInversionFixIterations) {
// If the matrix was not invertible, modify it a bit
for (int i = 0; i < n; i ++) {
coefficientsSave[i*n+i] = (1 + ((double)i)/((double)n)) * coefficientsSave[i*n+i];
}
int inverseSecondResult = Matrix::ArrayInverse(coefficientsSave, n, n);
assert(inverseSecondResult >= 0);
inverseResult = Matrix::ArrayInverse(coefficientsSave, n, n);
numberOfMatrixModifications++;
}
if (inverseResult < 0) {
return - 1;
}
if (numberOfMatrixModifications > 0) {
for (int i = 0; i < n*n; i++) {
coefficients[i] = coefficientsSave[i];
}
}
Multiplication::computeOnArrays<double>(coefficients, constants, solutions, n, n, 1);
return 0;
}
}

3
apps/regression/model/model.h
View File

@@ -42,6 +42,7 @@ private:
// Levenberg-Marquardt
static constexpr double k_maxIterations = 1000;
static constexpr double k_maxMatrixInversionFixIterations = 10;
static constexpr double k_initialLambda = 0.001;
static constexpr double k_lambdaFactor = 10;
static constexpr double k_chi2ChangeCondition = 0.001;
@@ -52,7 +53,7 @@ private:
double alphaPrimeCoefficient(Store * store, int series, double * modelCoefficients, int k, int l, double lambda) const;
double alphaCoefficient(Store * store, int series, double * modelCoefficients, int k, int l) const;
double betaCoefficient(Store * store, int series, double * modelCoefficients, int k) const;
void solveLinearSystem(double * solutions, double * coefficients, double * constants, int solutionDimension, Poincare::Context * context);
int solveLinearSystem(double * solutions, double * coefficients, double * constants, int solutionDimension, Poincare::Context * context);
};
}