Skip to content

Commit

Permalink
Factored LOO and K-Fold into MetaModelValidation
Browse files Browse the repository at this point in the history
  • Loading branch information
@mbaudin47
mbaudin47 committed May 5, 2024
1 parent 24e3b97 commit b4a6ce9
Show file tree
Hide file tree
Showing 4 changed files with 135 additions and 93 deletions.
55 changes: 8 additions & 47 deletions lib/src/Uncertainty/Algorithm/MetaModel/FunctionalChaos/FunctionalChaosValidation.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -116,18 +116,13 @@ SplitterImplementation FunctionalChaosValidation::getSplitter() const
/* Compute cross-validation Leave-One-Out metamodel predictions */
Sample FunctionalChaosValidation::ComputeMetamodelLeaveOneOutPredictions(
const FunctionalChaosResult & functionalChaosResult,
const LeaveOneOutSplitter &)
const LeaveOneOutSplitter & splitter)
{
const Sample outputSample(functionalChaosResult.getOutputSample());
const Sample residualsSample(functionalChaosResult.getSampleResiduals());
const Sample inputSample(functionalChaosResult.getInputSample());
const UnsignedInteger sampleSize = inputSample.getSize();
const FunctionCollection reducedBasis(functionalChaosResult.getReducedBasis());
const UnsignedInteger reducedBasisSize = reducedBasis.getSize();
if (reducedBasisSize > sampleSize)
throw InvalidArgumentException(HERE) << "Error: the sample size is: " << sampleSize <<
" which is lower than the basis size: " << reducedBasisSize;
const Sample outputSample(functionalChaosResult.getOutputSample());
const UnsignedInteger outputDimension = outputSample.getDimension();
const Function transformation(functionalChaosResult.getTransformation());
const Sample standardSample(transformation(inputSample));
DesignProxy designProxy(standardSample, reducedBasis);
Expand All @@ -137,17 +132,9 @@ Sample FunctionalChaosValidation::ComputeMetamodelLeaveOneOutPredictions(
const String methodName(ResourceMap::GetAsString("LeastSquaresExpansion-DecompositionMethod"));
LeastSquaresMethod leastSquaresMethod(LeastSquaresMethod::Build(methodName, designProxy, allIndices));
leastSquaresMethod.update(Indices(0), allIndices, Indices(0));
const Point diagonalH = leastSquaresMethod.getHDiag();
Sample cvPredictions(sampleSize, outputDimension);
for (UnsignedInteger j = 0; j < outputDimension; ++j)
for (UnsignedInteger i = 0; i < sampleSize; ++i)
{
if (diagonalH[i] == 1.0)
throw InvalidArgumentException(HERE)
<< "The leverage of observation #" << i
<< " is equal to 1. Cannot divide by zero.";
cvPredictions(i, j) = outputSample(i, j) - residualsSample(i, j) / (1.0 - diagonalH[i]);
} // For observations indices
const Point hMatrixDiag = leastSquaresMethod.getHDiag();
const Sample cvPredictions(MetaModelValidation::ComputeMetamodelLeaveOneOutPredictions(
outputSample, residualsSample, hMatrixDiag, splitter));
return cvPredictions;
}

Expand All @@ -156,15 +143,11 @@ Sample FunctionalChaosValidation::ComputeMetamodelKFoldPredictions(
const FunctionalChaosResult & functionalChaosResult,
const KFoldSplitter & splitter)
{
const Sample outputSample(functionalChaosResult.getOutputSample());
const Sample residualsSample(functionalChaosResult.getSampleResiduals());
const Sample inputSample(functionalChaosResult.getInputSample());
const UnsignedInteger sampleSize = inputSample.getSize();
const FunctionCollection reducedBasis(functionalChaosResult.getReducedBasis());
const UnsignedInteger reducedBasisSize = reducedBasis.getSize();
if (reducedBasisSize > sampleSize)
throw InvalidArgumentException(HERE) << "Error: the sample size is: " << sampleSize <<
" which is lower than the basis size: " << reducedBasisSize;
const Sample outputSample(functionalChaosResult.getOutputSample());
const UnsignedInteger outputDimension = outputSample.getDimension();
const Function transformation(functionalChaosResult.getTransformation());
const Sample standardSample(transformation(inputSample));
Expand All @@ -175,31 +158,9 @@ Sample FunctionalChaosValidation::ComputeMetamodelKFoldPredictions(
const String methodName(ResourceMap::GetAsString("LeastSquaresExpansion-DecompositionMethod"));
LeastSquaresMethod leastSquaresMethod(LeastSquaresMethod::Build(methodName, designProxy, allIndices));
leastSquaresMethod.update(Indices(0), allIndices, Indices(0));
Sample cvPredictions(sampleSize, outputDimension);
const SymmetricMatrix projectionMatrix(leastSquaresMethod.getH());
// Compute K-Fold error
UnsignedInteger kParameter = splitter.getSize();
Indices indicesTest;
for (UnsignedInteger foldIndex = 0; foldIndex < kParameter; ++foldIndex)
{
splitter.generate(indicesTest);
const UnsignedInteger foldSize = indicesTest.getSize();
SymmetricMatrix projectionKFoldMatrix(foldSize);
for (UnsignedInteger i1 = 0; i1 < foldSize; ++i1)
for (UnsignedInteger i2 = 0; i2 < 1 + i1; ++i2)
projectionKFoldMatrix(i1, i2) = projectionMatrix(indicesTest[i1], indicesTest[i2]);
const IdentityMatrix identityMatrix(foldSize);
const SymmetricMatrix reducedMatrix(identityMatrix - projectionKFoldMatrix);
const Sample residualsSampleKFoldTest(residualsSample.select(indicesTest));
Matrix multipleRightHandSide(foldSize, outputDimension);
for (UnsignedInteger j = 0; j < outputDimension; ++j)
for (UnsignedInteger i = 0; i < foldSize; ++i)
multipleRightHandSide(i, j) = residualsSampleKFoldTest(i, j);
const Matrix residualsKFoldMatrix(reducedMatrix.solveLinearSystem(multipleRightHandSide));
for (UnsignedInteger j = 0; j < outputDimension; ++j)
for (UnsignedInteger i = 0; i < foldSize; ++i)
cvPredictions(indicesTest[i], j) = outputSample(indicesTest[i], j) - residualsKFoldMatrix(i, j);
} // For fold indices
const Sample cvPredictions(MetaModelValidation::ComputeMetamodelKFoldPredictions(
outputSample, residualsSample, projectionMatrix, splitter));
return cvPredictions;
}

Expand Down
55 changes: 12 additions & 43 deletions lib/src/Uncertainty/Algorithm/MetaModel/LinearModel/LinearModelValidation.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -100,58 +100,27 @@ String LinearModelValidation::__repr__() const

/* Compute cross-validation predictions */
Sample LinearModelValidation::ComputeMetamodelLeaveOneOutPredictions(
const LinearModelResult & linearModelResult, const LeaveOneOutSplitter &)
const LinearModelResult & linearModelResult,
const LeaveOneOutSplitter & splitter)
{
// The residual is ri = g(xi) - tilde{g}(xi) where g is the model
// and tilde(g) is the metamodel.
// Hence the metamodel prediction is tilde{g}(xi) = yi - ri.
const Sample residualsSample(linearModelResult.getSampleResiduals());
const UnsignedInteger sampleSize = residualsSample.getSize();
const UnsignedInteger basisSize = linearModelResult.getBasis().getSize();
if (basisSize > sampleSize)
throw InvalidArgumentException(HERE) << "Error: the sample size is: " << sampleSize <<
" which is lower than the basis size: " << basisSize;
const Sample outputSample(linearModelResult.getOutputSample());
Sample cvPredictions(sampleSize, 1);
const Point diagonalH(linearModelResult.getLeverages());
for (UnsignedInteger i = 0; i < sampleSize; ++i)
cvPredictions(i, 0) = outputSample(i, 0) - residualsSample(i, 0) / (1.0 - diagonalH[i]);
const Sample residualsSample(linearModelResult.getSampleResiduals());
const Point hMatrixDiag(linearModelResult.getLeverages());
const Sample cvPredictions(MetaModelValidation::ComputeMetamodelLeaveOneOutPredictions(
outputSample, residualsSample, hMatrixDiag, splitter));
return cvPredictions;
}

/* Compute cross-validation predictions */
Sample LinearModelValidation::ComputeMetamodelKFoldPredictions(
const LinearModelResult & linearModelResult, const KFoldSplitter & splitter)
const LinearModelResult & linearModelResult,
const KFoldSplitter & splitter)
{
// The residual is ri = g(xi) - tilde{g}(xi) where g is the model
// and tilde(g) is the metamodel.
// Hence the metamodel prediction is tilde{g}(xi) = yi - ri.
const Sample residualsSample(linearModelResult.getSampleResiduals());
const UnsignedInteger sampleSize = residualsSample.getSize();
const UnsignedInteger basisSize = linearModelResult.getBasis().getSize();
if (basisSize > sampleSize)
throw InvalidArgumentException(HERE) << "Error: the sample size is: " << sampleSize <<
" which is lower than the basis size: " << basisSize;
const Sample outputSample(linearModelResult.getOutputSample());
Sample cvPredictions(sampleSize, 1);
LeastSquaresMethod lsMethod(linearModelResult.buildMethod());
SymmetricMatrix projectionMatrix(lsMethod.getH());
UnsignedInteger kParameter = splitter.getSize();
for (UnsignedInteger foldIndex = 0; foldIndex < kParameter; ++foldIndex)
{
Indices indicesTest;
splitter.generate(indicesTest);
const UnsignedInteger foldSize = indicesTest.getSize();
SymmetricMatrix projectionKFoldMatrix(foldSize);
for (UnsignedInteger i1 = 0; i1 < foldSize; ++i1)
for (UnsignedInteger i2 = 0; i2 < 1 + i1; ++i2)
projectionKFoldMatrix(i1, i2) = projectionMatrix(indicesTest[i1], indicesTest[i2]);
const SymmetricMatrix reducedMatrix(IdentityMatrix(foldSize) - projectionKFoldMatrix);
const Sample residualsSampleKFoldTest(residualsSample.select(indicesTest));
const Point residualsKFold(reducedMatrix.solveLinearSystem(residualsSampleKFoldTest.asPoint()));
for (UnsignedInteger i = 0; i < foldSize; ++i)
cvPredictions(indicesTest[i], 0) = outputSample(indicesTest[i], 0) - residualsKFold[i];
} // Loop over the folds
const Sample residualsSample(linearModelResult.getSampleResiduals());
SymmetricMatrix projectionMatrix(linearModelResult.buildMethod().getH());
const Sample cvPredictions(MetaModelValidation::ComputeMetamodelKFoldPredictions(
outputSample, residualsSample, projectionMatrix, splitter));
return cvPredictions;
}

Expand Down
101 changes: 101 additions & 0 deletions lib/src/Uncertainty/Algorithm/MetaModel/MetaModelValidation.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -173,6 +173,107 @@ GridLayout MetaModelValidation::drawValidation() const
return grid;
}

/* Compute cross-validation predictions */
Sample MetaModelValidation::ComputeMetamodelLeaveOneOutPredictions(
const Sample & outputSample,
const Sample & residual,
const Point & hMatrixDiag,
const LeaveOneOutSplitter & splitter)
{
// The residual is ri = g(xi) - tilde{g}(xi) where g is the model
// and tilde(g) is the metamodel.
// Hence the metamodel prediction is tilde{g}(xi) = yi - ri.
const UnsignedInteger sampleSize = outputSample.getSize();
const UnsignedInteger outputDimension = outputSample.getDimension();
if (residual.getSize() != sampleSize)
throw InvalidArgumentException(HERE)
<< "Error: the residual sample size is: "
<< residual.getSize() <<
" but the output sample size is: " << sampleSize;
if (residual.getDimension() != outputSample.getDimension())
throw InvalidArgumentException(HERE)
<< "Error: the output sample dimension is: "
<< outputSample.getDimension()
<< " which is different from the residual output dimension: "
<< residual.getDimension();
if (splitter.getN() != sampleSize)
throw InvalidArgumentException(HERE)
<< "Error: the splitter size is: " << splitter.getSize() <<
" but the output sample size is " << sampleSize;
if (hMatrixDiag.getDimension() != sampleSize)
throw InvalidArgumentException(HERE)
<< "Error: the H matrix diagional dimension is: " << hMatrixDiag.getDimension() <<
" but the output sample size is " << sampleSize;
Sample cvPredictions(sampleSize, outputDimension);
for (UnsignedInteger j = 0; j < outputDimension; ++j)
for (UnsignedInteger i = 0; i < sampleSize; ++i)
{
if (hMatrixDiag[i] == 1.0)
throw InvalidArgumentException(HERE)
<< "The leverage of observation #" << i
<< " is equal to 1. Cannot divide by zero.";
cvPredictions(i, j) = outputSample(i, j) - residual(i, j) / (1.0 - hMatrixDiag[i]);
} // For observations indices
return cvPredictions;
}

/* Compute cross-validation predictions */
Sample MetaModelValidation::ComputeMetamodelKFoldPredictions(
const Sample & outputSample,
const Sample & residual,
const SymmetricMatrix & projectionMatrix,
const KFoldSplitter & splitter)
{
// The residual is ri = g(xi) - tilde{g}(xi) where g is the model
// and tilde(g) is the metamodel.
// Hence the metamodel prediction is tilde{g}(xi) = yi - ri.
const UnsignedInteger sampleSize = outputSample.getSize();
const UnsignedInteger outputDimension = outputSample.getDimension();
if (residual.getSize() != sampleSize)
throw InvalidArgumentException(HERE)
<< "Error: the residual sample size is: "
<< residual.getDimension() <<
" but the output sample size is: " << sampleSize;
if (residual.getDimension() != outputSample.getDimension())
throw InvalidArgumentException(HERE)
<< "Error: the output sample dimension is: "
<< outputSample.getDimension()
<< " which is different from the residual output dimension: "
<< residual.getDimension();
if (splitter.getN() != sampleSize)
throw InvalidArgumentException(HERE)
<< "Error: the splitter size is: " << splitter.getSize() <<
" but the output sample size is " << sampleSize;
if (projectionMatrix.getDimension() != sampleSize)
throw InvalidArgumentException(HERE)
<< "Error: the H matrix diagional dimension is: " << projectionMatrix.getDimension() <<
" but the output sample size is " << sampleSize;
Sample cvPredictions(sampleSize, outputDimension);
UnsignedInteger kParameter = splitter.getSize();
Indices indicesTest;
for (UnsignedInteger foldIndex = 0; foldIndex < kParameter; ++foldIndex)
{
splitter.generate(indicesTest);
const UnsignedInteger foldSize = indicesTest.getSize();
SymmetricMatrix projectionKFoldMatrix(foldSize);
for (UnsignedInteger i1 = 0; i1 < foldSize; ++i1)
for (UnsignedInteger i2 = 0; i2 < 1 + i1; ++i2)
projectionKFoldMatrix(i1, i2) = projectionMatrix(indicesTest[i1], indicesTest[i2]);
const IdentityMatrix identityMatrix(foldSize);
const SymmetricMatrix reducedMatrix(identityMatrix - projectionKFoldMatrix);
const Sample residualsSampleKFoldTest(residual.select(indicesTest));
Matrix multipleRightHandSide(foldSize, outputDimension);
for (UnsignedInteger j = 0; j < outputDimension; ++j)
for (UnsignedInteger i = 0; i < foldSize; ++i)
multipleRightHandSide(i, j) = residualsSampleKFoldTest(i, j);
const Matrix residualsKFoldMatrix(reducedMatrix.solveLinearSystem(multipleRightHandSide));
for (UnsignedInteger j = 0; j < outputDimension; ++j)
for (UnsignedInteger i = 0; i < foldSize; ++i)
cvPredictions(indicesTest[i], j) = outputSample(indicesTest[i], j) - residualsKFoldMatrix(i, j);
} // For fold indices
return cvPredictions;
}

/* Method save() stores the object through the StorageManager */
void MetaModelValidation::save(Advocate & adv) const
{
Expand Down
17 changes: 14 additions & 3 deletions lib/src/Uncertainty/Algorithm/MetaModel/openturns/MetaModelValidation.hxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,11 +24,11 @@
#include "openturns/Function.hxx"
#include "openturns/Distribution.hxx"
#include "openturns/GridLayout.hxx"
#include "openturns/LeaveOneOutSplitter.hxx"
#include "openturns/KFoldSplitter.hxx"

BEGIN_NAMESPACE_OPENTURNS



/**
* @class MetaModelValidation
*
Expand Down Expand Up @@ -98,8 +98,19 @@ protected:
/** R2 score */
mutable Point r2Score_;

private:
/** Compute cross-validation leave-one-out predictions */
static Sample ComputeMetamodelLeaveOneOutPredictions(const Sample & outputSample,
const Sample & residual,
const Point & leverages,
const LeaveOneOutSplitter & splitter);

/** Compute cross-validation K-Fold predictions */
static Sample ComputeMetamodelKFoldPredictions(const Sample & outputSample,
const Sample & residual,
const SymmetricMatrix & projectionMatrix,
const KFoldSplitter & splitter);

private:

}; /* class MetaModelValidation */

Expand Down

0 comments on commit b4a6ce9

Please sign in to comment.