Changeset cec7a5

Timestamp:

Oct 6, 2011, 7:17:47 AM (14 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

59e063

Parents:

7a0340

git-author:

Frederik Heber <heber@…> (04/23/11 20:20:38)

git-committer:

Frederik Heber <heber@…> (10/06/11 07:17:47)

Message:

Added Singular Value Decomposition to MatrixContent.

• this allows for solving over-determined linear system of equations in a least

squares sense:

Ax - b

2.

• Unit test added on SVD and solving of an example system.

Location:

LinearAlgebra/src

Files:

• LinearAlgebra/MatrixContent.cpp (modified) (2 diffs)
• LinearAlgebra/MatrixContent.hpp (modified) (1 diff)
• unittests/MatrixContentUnitTest.cpp (modified) (3 diffs)
• unittests/MatrixContentUnitTest.hpp (modified) (2 diffs)

LinearAlgebra/src/LinearAlgebra/MatrixContent.cpp

@@ -86 +86 @@
 }
 
-/** Constructor that parses square matrix from a stream.
+/** Constructor that parses matrix from a stream.
  *
  * \note Matrix dimensions can be preparsed via
  * MatrixContent::preparseMatrixDimensions() without harming the stream.
  *
+ * \param _row number of rows
+ * \param _column number of columns
  * \param &inputstream stream to parse from
  */
@@ -509 +511 @@
     }
   return *this;
+}
+
+/** Performs a SVD on the contained matrix.
+ *
+ * A SVD is the factorization of a MxN matrix A into the product form \f$U S V^t\f$, where
+ *  -# U is a MxN orthogonal matrix
+ *  -# S is a NxN diagonal matrix (here just a vector)
+ *  -# V is a NxN orthogonal, square matrix.
+ *
+ * The SVD is only done if MatrixContent::rows is larger equal to MatrixContent::columns, i.e.
+ * N >= M.
+ *
+ * \note The matrix content of this matrix, \a V and \a S are overwritten in the process.
+ *
+ * \param V orthogonal square matrix on return
+ * \param S diagonal matrix on return
+ */
+void MatrixContent::performSingularValueDecomposition(MatrixContent &V, VectorContent &S)
+{
+  ASSERT(rows >= columns,
+      "MatrixContent::performSingularValueDecomposition() - row count "
+      +toString(rows)+" must be larger than or equal to column count "+toString(columns)+".");
+  ASSERT(V.getRows() == V.getColumns(),
+      "MatrixContent::performSingularValueDecomposition() - Matrix V must be square!");
+  ASSERT(V.getRows() == columns,
+      "MatrixContent::performSingularValueDecomposition() - Given matrix V does not have the right dimension "
+      +toString(V.getRows())+"!="+toString(rows)+".");
+  ASSERT(S.getDimension() == columns,
+      "MatrixContent::performSingularValueDecomposition() - Vector S does not haver the right dimension "
+      +toString(S.getDimension())+"!="+toString(columns)+".");
+  gsl_matrix *A = content;
+  gsl_vector *work = gsl_vector_alloc(columns);
+  gsl_linalg_SV_decomp(A,V.content,S.content,work);
+  gsl_vector_free(work);
+}
+
+/** Performs a SVD on the contained matrix and solves the linear system \f$Ax=b\f$ with given \a b.
+ *
+ * The SVD is only done if MatrixContent::rows is larger equal to MatrixContent::columns.
+ *
+ * \note The matrix content is overwritten in the process.
+ *
+ * \param V NxN orthogonal, diagonal matrix from SVD decomposition
+ * \param S N vector representing the diagonal elements from SVD decomposition
+ * \param *b right-hand side of linear equation \f$Ax=b\f$ to solve
+ * \return *x solution, if rows>columns this will be the solution in the least squares sense \f$||Ax-b||^2\f$
+ */
+VectorContent MatrixContent::solveOverdeterminedLinearEquation(MatrixContent &V, VectorContent &S, const VectorContent &b)
+{
+  ASSERT(b.getDimension() == rows,
+      "MatrixContent::performSingularValueDecomposition() - Given vector b does not have the right dimension "
+      +toString(b.getDimension())+"!="+toString(rows)+".");
+  ASSERT(V.getRows() == V.getColumns(),
+      "MatrixContent::performSingularValueDecomposition() - Matrix V must be square!");
+  ASSERT(V.getRows() == columns,
+      "MatrixContent::performSingularValueDecomposition() - Given matrix V does not have the right dimension "
+      +toString(V.getRows())+"!="+toString(rows)+".");
+  ASSERT(S.getDimension() == columns,
+      "MatrixContent::performSingularValueDecomposition() - Vector S does not haver the right dimension "
+      +toString(S.getDimension())+"!="+toString(columns)+".");
+
+  gsl_vector *x = gsl_vector_alloc(columns);
+  gsl_linalg_SV_solve(content,V.content,S.content,b.content,x);
+  return VectorContent(x, true);
+}
+
+/** Performs a SVD on the contained matrix and solves the linear system \f$Ax=b\f$ with given \a b.
+ *
+ * The SVD is only done if MatrixContent::rows is larger equal to MatrixContent::columns.
+ *
+ * \note The matrix content is overwritten in the process.
+ *
+ * \param *b right-hand side of linear equation \f$Ax=b\f$ to solve.
+ * \return *x solution, if rows>columns this will be the solution in the least squares sense \f$||Ax-b||^2\f$
+ */
+VectorContent MatrixContent::solveOverdeterminedLinearEquation(const VectorContent &b)
+{
+  ASSERT(b.getDimension() == rows,
+      "MatrixContent::performSingularValueDecomposition() - Given vector b does not have the right dimension "
+      +toString(b.getDimension())+"!="+toString(rows)+".");
+  MatrixContent V(columns, columns);
+  VectorContent S(columns);
+  performSingularValueDecomposition(V,S);
+  return solveOverdeterminedLinearEquation(V,S,b);
 }
 

LinearAlgebra/src/LinearAlgebra/MatrixContent.hpp

@@ -98 +98 @@
   gsl_vector* transformToEigenbasis();
   void sortEigenbasis(gsl_vector *eigenvalues);
+  void performSingularValueDecomposition(MatrixContent &V, VectorContent &S);
+  VectorContent solveOverdeterminedLinearEquation(const VectorContent &b);
+  VectorContent solveOverdeterminedLinearEquation(MatrixContent &V, VectorContent &S, const VectorContent &b);
 
   // Properties

LinearAlgebra/src/unittests/MatrixContentUnitTest.cpp

@@ -24 +24 @@
 #include <cppunit/ui/text/TestRunner.h>
 
+#include <cmath>
+#include <limits>
+
 #include "MatrixContentUnitTest.hpp"
 
 #include "MatrixContent.hpp"
+#include "VectorContent.hpp"
 
 #ifdef HAVE_TESTRUNNER
@@ -33 +37 @@
 
 /********************************************** Test classes **************************************/
+
+const std::string matrixA = "   9.962848439806617e-01   5.950413028139907e-01   3.172650036465889e-01\n\
+   2.070939109924951e-01   8.365712117773689e-01   3.980567886073226e-01\n\
+   1.021600573071414e-01   8.061359922326598e-02   8.493537275043391e-01\n\
+   6.955361514014282e-01   6.873206387346102e-02   2.649165458481005e-01";
+const std::string vectorS = "1.630290761001098e+00                    7.624807744374841e-01                    6.374093742147465e-01";
 
 // Registers the fixture into the 'registry'
@@ -293 +303 @@
   delete n;
 }
+
+/** Unit test for MatrixContent::performSingularValueDecomposition().
+ *
+ */
+void MatrixContentTest::SVDTest()
+{
+  // setup "A", U,S,V
+  std::stringstream Astream(matrixA);
+  std::stringstream Sstream(vectorS);
+  MatrixContent A(m->getRows(), m->getColumns(), Astream);
+  VectorContent S_expected(m->getColumns(), Sstream);
+  *m = A;
+
+  // check SVD
+  VectorContent S(m->getColumns());
+  MatrixContent V(m->getColumns(), m->getColumns());
+  A.performSingularValueDecomposition(V,S);
+  MatrixContent S_diag(m->getColumns(), m->getColumns());
+  for (size_t index = 0; index < m->getColumns(); ++index)
+    S_diag.set(index, index, S[index]);
+  CPPUNIT_ASSERT_EQUAL(*m,A*S_diag*(V.transpose()));
+  for (size_t index = 0; index < m->getColumns(); ++index) {
+    //CPPUNIT_ASSERT_EQUAL(S_expected[index], S[index]);
+    CPPUNIT_ASSERT(fabs(S_expected[index] - S[index]) < numeric_limits<double>::epsilon()*10.);
+  }
+
+  // set up right-hand side
+  VectorContent b(4);
+  b[0] = 1.;
+  b[1] = 0.;
+  b[2] = 0.;
+  b[3] = 0.;
+
+  // SVD
+  VectorContent x_expected(3);
+  x_expected[0] = -0.00209169;
+  x_expected[1] = -0.325399;
+  x_expected[2] = 0.628004;
+  const VectorContent x_result(A.solveOverdeterminedLinearEquation(V,S,b));
+
+  // check result
+  for (size_t index = 0; index < m->getColumns(); ++index) {
+    CPPUNIT_ASSERT(fabs(x_expected[index] - x_result[index]) < 1.e-6);
+  }
+}

LinearAlgebra/src/unittests/MatrixContentUnitTest.hpp

@@ -30 +30 @@
     CPPUNIT_TEST (PropertiesTest );
     CPPUNIT_TEST (ReadWriteTest );
+    CPPUNIT_TEST (SVDTest );
     CPPUNIT_TEST_SUITE_END();
 
@@ -42 +43 @@
     void PropertiesTest();
     void ReadWriteTest();
+    void SVDTest();
 
 private: