From c56652cda9bc58959ddfb09ac839f03f434c0bf7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9a=20Saviot?= Date: Fri, 15 Jun 2018 17:27:07 +0200 Subject: [PATCH] [apps/reg] Algorithm more robust against non invertible matrices --- apps/regression/model/model.cpp | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/apps/regression/model/model.cpp b/apps/regression/model/model.cpp index a002728d4..def197f1a 100644 --- a/apps/regression/model/model.cpp +++ b/apps/regression/model/model.cpp @@ -148,8 +148,23 @@ double Model::betaCoefficient(Store * store, int series, double * modelCoefficie void 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]; + for (int i = 0; i < n * n; i++) { + coefficientsSave[i] = coefficients[i]; + } int inverseResult = Matrix::ArrayInverse(coefficients, n, n); - assert(inverseResult >= 0); + if (inverseResult < 0) { + // 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); + for (int i = 0; i < n*n; i++) { + coefficients[i] = coefficientsSave[i]; + } + } Multiplication::computeOnArrays(coefficients, constants, solutions, n, n, 1); }