From e13ef0a41f45cc7e6de5db39262384f825a7c169 Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Thu, 12 Apr 2018 15:38:42 +0200 Subject: [PATCH] Fix issue with AD SymmetricTensor eigenvalue/vector calculations. For auto-differentiable numbers, when the input tensor was diagonal the sensitivities of the eigenvalues/vector with respect to one another were deduced. This means that although the correct results were being returned, their derivatives were incorrect. This patch attempts to correct this by perturbing a diagonal input tensor (of AD number type) and computing an approximation to the eigenvalues/vectors. --- .../changes/minor/20180412Jean-PaulPelteret-3 | 9 + include/deal.II/base/symmetric_tensor.h | 43 +-- .../deal.II/base/symmetric_tensor.templates.h | 262 +++++++++++++++++- 3 files changed, 267 insertions(+), 47 deletions(-) create mode 100644 doc/news/changes/minor/20180412Jean-PaulPelteret-3 diff --git a/doc/news/changes/minor/20180412Jean-PaulPelteret-3 b/doc/news/changes/minor/20180412Jean-PaulPelteret-3 new file mode 100644 index 0000000000..b26c208a4d --- /dev/null +++ b/doc/news/changes/minor/20180412Jean-PaulPelteret-3 @@ -0,0 +1,9 @@ +Fixed: When using diagonal SymmetricTensors with automatic-differentiable numbers, computations +using the eigenvalue() and eigenvector() functions would return correct values of the +eigenvalues/vectors. However, the derivatives of these values/vectors were incorrect as +the sensitivities of the eigenvalues/vectors with respect to one another was not encoded +in the returned result. Therefore, under these specific conditions the returned result is now +a more coarse approximation for the eigenvalues/vectors but with the trade-off that a +meaningful approximation of the derivative of the result can now be computed. +
+(Jean-Paul Pelteret, 2018/04/12) diff --git a/include/deal.II/base/symmetric_tensor.h b/include/deal.II/base/symmetric_tensor.h index baf634f95e..08b857a1d2 100644 --- a/include/deal.II/base/symmetric_tensor.h +++ b/include/deal.II/base/symmetric_tensor.h @@ -1313,7 +1313,7 @@ namespace internal Number diagonal_sum = internal::NumberType::value(0.0); for (unsigned int i=0; i(N); (void)typical_diagonal_element; unsigned int p[N]; @@ -3112,22 +3112,6 @@ enum struct SymmetricTensorEigenvectorMethod -/** - * Return the eigenvalues and eigenvectors of a symmetric 1x1 tensor of rank 2. - * - * @relatesalso SymmetricTensor - * @author Jean-Paul Pelteret, 2017 - */ -template -std::array >,1> -eigenvectors (const SymmetricTensor<2,1,Number> &T, - const SymmetricTensorEigenvectorMethod /*method*/ = SymmetricTensorEigenvectorMethod::ql_implicit_shifts) -{ - return { {std::make_pair(T[0][0], Tensor<1,1,Number>({1.0}))} }; -} - - - /** * Return the eigenvalues and eigenvectors of a real-valued rank-2 symmetric * tensor $T$. The array of matched eigenvalue and eigenvector pairs is sorted @@ -3158,30 +3142,7 @@ eigenvectors (const SymmetricTensor<2,1,Number> &T, template std::array >,dim> eigenvectors (const SymmetricTensor<2,dim,Number> &T, - const SymmetricTensorEigenvectorMethod method = SymmetricTensorEigenvectorMethod::ql_implicit_shifts) -{ - std::array >,dim> eig_vals_vecs; - - switch (method) - { - case SymmetricTensorEigenvectorMethod::hybrid: - eig_vals_vecs = internal::SymmetricTensorImplementation::hybrid(T); - break; - case SymmetricTensorEigenvectorMethod::ql_implicit_shifts: - eig_vals_vecs = internal::SymmetricTensorImplementation::ql_implicit_shifts(T); - break; - case SymmetricTensorEigenvectorMethod::jacobi: - eig_vals_vecs = internal::SymmetricTensorImplementation::jacobi(T); - break; - default: - AssertThrow(false, ExcNotImplemented()); - } - - // Sort in descending order before output. - std::sort(eig_vals_vecs.begin(), eig_vals_vecs.end(), - internal::SymmetricTensorImplementation::SortEigenValuesVectors()); - return eig_vals_vecs; -} + const SymmetricTensorEigenvectorMethod method = SymmetricTensorEigenvectorMethod::ql_implicit_shifts); diff --git a/include/deal.II/base/symmetric_tensor.templates.h b/include/deal.II/base/symmetric_tensor.templates.h index e8f80e8798..0f5cc7c282 100644 --- a/include/deal.II/base/symmetric_tensor.templates.h +++ b/include/deal.II/base/symmetric_tensor.templates.h @@ -17,7 +17,10 @@ #define dealii_symmetric_tensor_templates_h +#include #include +#include +#include #include @@ -62,8 +65,8 @@ eigenvalues (const SymmetricTensor<2,2,Number> &T) const std::array eig_vals = { { - static_cast(0.5*(tr_T + sqrt_desc)), - static_cast(0.5*(tr_T - sqrt_desc)) + internal::NumberType::value(0.5*(tr_T + sqrt_desc)), + internal::NumberType::value(0.5*(tr_T - sqrt_desc)) } }; Assert(eig_vals[0] >= eig_vals[1], ExcMessage("The eigenvalue ordering is incorrect.")); @@ -110,11 +113,14 @@ eigenvalues (const SymmetricTensor<2,3,Number> &T) // The value of tmp_2 should be within [-1,1], however // floating point errors might place it slightly outside - // this range. It is therefore necessary to correct for it + // this range. It is therefore necessary to correct for it. + // Note: The three results in the conditional may lead to different + // number types when using Sacado numbers, so we cast them when + // necessary to a consistent result type. const Number phi = - (tmp_2 <= -1.0 ? Number(numbers::PI/3.0) : - (tmp_2 >= 1.0 ? Number(0.0) : - std::acos(tmp_2)/3.0)); + (tmp_2 <= -1.0 ? internal::NumberType::value(numbers::PI/3.0) : + (tmp_2 >= 1.0 ? internal::NumberType::value(0.0) : + internal::NumberType::value(std::acos(tmp_2)/3.0))); // Due to the trigonometric solution, the computed eigenvalues // should be predictably in the order eig1 >= eig2 >= eig3... @@ -704,11 +710,255 @@ namespace internal return eig_vals_vecs; } + + + template + Tensor<2,1,Number> + dediagonalize_tensor (const dealii::SymmetricTensor<2,1,Number> &T, + const double /*rotation_angle*/, + const unsigned int /*axis*/ = 0) + { + AssertThrow(false, ExcNotImplemented()); + return Tensor<2,1,Number> ({{T[0][0]}}); + } + + + template + Tensor<2,2,Number> + dediagonalize_tensor (const dealii::SymmetricTensor<2,2,Number> &T, + const double rotation_angle, + const unsigned int /*axis*/ = 0) + { + const Tensor<2,2> R = dealii::Physics::Transformations::Rotations::rotation_matrix_2d(rotation_angle); + return R*T; + } + + + template + Tensor<2,3,Number> + dediagonalize_tensor (const dealii::SymmetricTensor<2,3,Number> &T, + const double rotation_angle, + const unsigned int axis = 0) + { + Assert(axis<3, ExcIndexRange(axis,0,3)); + + Tensor<2,3> R; + switch (axis) + { + case (0): + R = dealii::Physics::Transformations::Rotations::rotation_matrix_3d({1,0,0},rotation_angle); + break; + case (1): + R = dealii::Physics::Transformations::Rotations::rotation_matrix_3d({0,1,0},rotation_angle); + break; + case (2): + R = dealii::Physics::Transformations::Rotations::rotation_matrix_3d({0,0,1},rotation_angle); + break; + default: + AssertThrow(false, ExcNotImplemented()); + break; + } + return R*T; + } + + + template + std::array >,1> + perform_eigenvector_decomposition (const SymmetricTensor<2,1,Number> &T, + const SymmetricTensorEigenvectorMethod /*method*/) + { + return { {std::make_pair(T[0][0], Tensor<1,1,Number>({1.0}))} }; + } + + + template + std::array >,dim> + perform_eigenvector_decomposition (const SymmetricTensor<2,dim,Number> &T, + const SymmetricTensorEigenvectorMethod method) + { + switch (method) + { + case SymmetricTensorEigenvectorMethod::hybrid: + return internal::SymmetricTensorImplementation::hybrid(T); + break; + case SymmetricTensorEigenvectorMethod::ql_implicit_shifts: + return internal::SymmetricTensorImplementation::ql_implicit_shifts(T); + break; + case SymmetricTensorEigenvectorMethod::jacobi: + return internal::SymmetricTensorImplementation::jacobi(T); + break; + default: + break; + } + + AssertThrow(false, ExcNotImplemented()); + return std::array >,dim> (); + } + + } // namespace SymmetricTensor } // namespace internal +template +std::array >,dim> +eigenvectors (const SymmetricTensor<2,dim,Number> &T, + const SymmetricTensorEigenvectorMethod method) +{ + // Not much to do when there's only a single entry + if (dim == 1) + return internal::SymmetricTensorImplementation::perform_eigenvector_decomposition(T,method); + + std::array >,dim> eig_vals_vecs; + + if (Differentiation::AD::is_ad_number::value && dim>1) + { + // If the tensor is diagonal, then we have a bit on an issue when using + // auto-differentiable numbers. The reason for this is that all of the + // algorithms have shortcuts by which to return result in this case. + // This artificially decouples the eigenvalues/vectors, which upon + // differentiation leads to the wrong result (each is, incorrectly, + // insensitive with respect to the other). To work around this manipulate + // tensor @p T in an objective manner: through an infinitesimal rotation we + // make it non-diagonal (although we introduce some numerical error). + bool is_diagonal = true; + for (unsigned int i=0; i::value) + { + // Adol-C taped + if (method == SymmetricTensorEigenvectorMethod::ql_implicit_shifts) + sf = (dim == 2 ? 2e11 : 2e11); + else if (method == SymmetricTensorEigenvectorMethod::jacobi) + sf = (dim == 2 ? 1e6 : 1e9); + else + AssertThrow(false, ExcNotImplemented()); + } + else if (Differentiation::AD::is_sacado_rad_number::value) + { + // Sacado::Rad + if (method == SymmetricTensorEigenvectorMethod::ql_implicit_shifts) + sf = (dim == 2 ? 1e8 : 1e9); + else if (method == SymmetricTensorEigenvectorMethod::jacobi) + sf = (dim == 2 ? 1e8 : 1e9); + else + AssertThrow(false, ExcNotImplemented()); + } + else + { + // Everything else + Assert(Differentiation::AD::is_tapeless_ad_number::value, ExcInternalError()); + Assert(Differentiation::AD::is_sacado_dfad_number::value || + Differentiation::AD::is_adolc_tapeless_number::value, ExcInternalError()); + + if (method == SymmetricTensorEigenvectorMethod::ql_implicit_shifts) + sf = (dim == 2 ? 1e7 : 2.5e7); + else if (method == SymmetricTensorEigenvectorMethod::jacobi) + sf = (dim == 2 ? 1e2 : 1e7); + else + AssertThrow(false, ExcNotImplemented()); + } + + typedef typename Differentiation::AD::ADNumberTraits::scalar_type scalar_type; + const double delta = sf*std::numeric_limits::epsilon(); + const double rotation_angle = delta*numbers::PI/180.0; + + Tensor<2,dim,Number> T_prime_ns; + if (dim == 2) + { + const Tensor<2,dim,Number> T_prime_ns = internal::SymmetricTensorImplementation::dediagonalize_tensor(T,rotation_angle); + + // We can't symmetrize the tensor, otherwise the sensitivities + // cancel out. So we take the upper triangle as an approximation + // instead. + // TODO[JPP]: Perform the eigen-decomposition on the non-symmetric T_prime_ns. + // This is, however, nontrivial to implement in this context. See: + // http://www.alglib.net/eigen/nonsymmetric/nonsymmetricevd.php + // https://groups.google.com/forum/#!topic/stan-users/QJe1TNioiyg + SymmetricTensor<2,dim,Number> T_prime; + for (unsigned int i=0; i T_prime; + Tensor<2,dim,Number> T_prime_ns; + for (unsigned int i=0; i()); + return eig_vals_vecs; +} + + + DEAL_II_NAMESPACE_CLOSE #endif -- 2.39.5