search

Open-Systems-Pharmacology / OSPSuite.Core

Core functionalities of the Open Systems Pharmacology Suite
Other
5 stars 8 forks source link

Improvements for calculateCovarianceMatrix #252

Open hdiedam opened 7 years ago

hdiedam commented 7 years ago

The code in calculateCovarianceMatrix can be improved.

In line with the original youtrack 47-7683: Although, the code used in this function currently is theoretically correct (because pinv(J'J)=inv(J'J)), it is numerically inferior and much more complex compared to the original suggestion. Using directly pinv(J)*S*pinv(J)', scalar times matrix multiplication instead of building the dense diagonal matrix and the API from https://numerics.mathdotnet.com/api/MathNet.Numerics.LinearAlgebra.Double/Matrix.htm should lead to (untestet):

private Matrix calculateCovarianceMatrix(JacobianMatrix jacobianMatrix, ResidualsResult residualsResult) { try { var pinv = jacobianMatrix.PseudoInverse(); var numberOfIdentifiedParameters = jacobianMatrix.ColumnCount; var numberOfResiduls = jacobianMatrix.RowCount; var s = residualsResult.SumResidual2 / (numberOfResiduls - numberOfIdentifiedParameters);

        return pinv.Multiply(s).TransposeAndMultiply(pinv);
     }
     catch (Exception e)
     {
        throw new MatrixCalculationException(Error.CovarianceMatrixCannotBeCalculated, e);
     }
  }

which should be a significant stability,readability and speed improvement.

@Yuri05 @KatrinCoboeken

msevestre commented 7 years ago

We have tests covering this function throughout the test suite. We could just swap the implementation and see if they are still working. I don't think speed improvement will be significant here since we do not calculate the covariate matrix very often. Improvement and stability however is well worth the change!

msevestre commented 6 years ago

@hdiedam JacobianMatrix is an internal object of OSPSuite and does not have a method PseudoInverse().

Could you please check the code and let me know what should be implemented? Going to remove it from 7.3.0 for now

msevestre commented 6 years ago

@hdiedam I guess I can simply create a Matrix object from .NET numeric filled with our matrix and go from there? It this what is meant?

hdiedam commented 6 years ago

Yes, I just assumed that JacobianMatrix was derived from Matrix, because the function returns a Matrix object and in line 81 you use "var fisherInv = fisher.PseudoInverse();". If JacobianMatrix itself is not a Matrix just create one using the underlying data array.