Pseudo-inverse of symmetric matrices using SVD / Utility for moving least squares #950

mrlag31 · 2023-09-11T20:19:40Z

In an attempt to break-up #946, this PR focuses only on the pseudo inverse used in the moving least squares algorithm. It includes the pseudo inverse itself in src/interpolation folder, tests and ~~an iterator for multi-dimensional views (for use in boost tests)~~ a macro to test equality on md views

dalg24 · 2023-09-12T11:55:11Z

src/interpolation/details/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp

+symmetricPseudoInverseSVDSerialKernel(InOutMatrix &io, SMatrix &s, UMatrix &u)
+{
+  using value_t = typename InOutMatrix::non_const_value_type;
+  std::size_t const size = io.extent(0);


Is io.extent(1) == size a precondition?
Presumably s and u must also be properly sized.

Yes and yes. io, s and u must have the same size and this should be guaranteed by the host function. But due to the fact I cannot use assert in kernels, I skipped checking it. (although, speaking of it, I could use KOKKOS_ASSERT)

src/interpolation/details/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp

aprokop · 2023-09-11T22:29:48Z

test/CMakeLists.txt

@@ -235,6 +235,12 @@ target_link_libraries(ArborX_Test_BoostAdapters.exe PRIVATE ArborX Boost::unit_t
 target_compile_definitions(ArborX_Test_BoostAdapters.exe PRIVATE BOOST_TEST_DYN_LINK)
 add_test(NAME ArborX_Test_BoostAdapters COMMAND ArborX_Test_BoostAdapters.exe)

+add_executable(ArborX_Test_InterpDetailsSymmPInvSVD.exe tstInterpDetailsSymmPInvSVD.cpp utf_main.cpp)


Probably change all the relevant names ArborX_Test_DetailsInterpolationSVD. This would emphasize that it is a Details test, as it starts with ArborX_Test_Details. We don't need to emphasize what exact SVD is done, so may omit that. And imho, Interpolation is better than Interp.

I will stick with InterpDetailsSVD because it mimics the location and name of the tested file (in src/interpolation/details). But maybe would it be better to have ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp elsewhere?

src/interpolation/details/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp

aprokop · 2023-09-13T15:52:48Z

src/interpolation/details/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp

+    value_t const a = sigma(p, p);
+    value_t const b = sigma(p, q);
+    value_t const c = sigma(q, q);


Suggested change

value_t const a = sigma(p, p);

value_t const b = sigma(p, q);

value_t const c = sigma(q, q);

auto const a = sigma(p, p);

auto const b = sigma(p, q);

auto const c = sigma(q, q);

aprokop · 2023-09-13T15:53:29Z

src/interpolation/details/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp

+    // | b | c |              | 0 | y |
+    // +---+---+              +---+---+
+
+    value_t cos, sin, x, y;


Probably a bad idea to name variables cos and sin. If one puts using namespace std; somewhere above, this will break.

aprokop · 2023-09-13T15:54:30Z

src/interpolation/details/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp

+    if (a == c)
+    {
+      cos = Kokkos::sqrt(value_t(2)) / 2;
+      sin = Kokkos::sqrt(value_t(2)) / 2;


Dunno if the compiler is smart enough to not do this computation twice.

Suggested change

sin = Kokkos::sqrt(value_t(2)) / 2;

sin = cos;

aprokop · 2023-09-13T15:55:09Z

src/interpolation/details/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp

+      value_t const u = (2 * b) / (a - c);
+      value_t const v = 1 / Kokkos::sqrt(u * u + 1);


Suggested change

value_t const u = (2 * b) / (a - c);

value_t const v = 1 / Kokkos::sqrt(u * u + 1);

auto const u = (2 * b) / (a - c);

auto const v = 1 / Kokkos::sqrt(u * u + 1);

aprokop · 2023-09-13T16:02:08Z

src/interpolation/details/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp

+
+template <typename ExecutionSpace, typename InvMatrices>
+void symmetricPseudoInverseSVD(ExecutionSpace const &space,
+                               InvMatrices &inv_matrices)


I think the name InvMatrices is misleading, given that it is both input and output.

I have renamed the type InOutMatrices, as the inverse is already implied by the function name.

src/interpolation/details/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp

dalg24 · 2023-09-14T17:31:38Z

src/interpolation/details/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp

+{
+  isSquareMatrix(mat);
+  using value_t = typename Matrix::non_const_value_type;
+  int const size = mat.extent(0);


Move this closer to the place you first use it.

dalg24 · 2023-09-14T22:53:13Z

src/interpolation/details/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp

+  struct
+  {
+    value_t max = 0;
+    int row = 0;
+    int col = 0;
+  } result;


I didn't even know you could return an object from a local unnamed class.
Sure why not given that these are implementation details.

test/ArborX_EnableViewComparison.hpp

masterleinad · 2023-09-18T19:20:05Z

src/interpolation/details/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp

+    int i = 0;
+    for (; i < p; i++)
+    {


Suggested change

int i = 0;

for (; i < p; i++)

{

for (int i = 0; i < p; i++)

{

masterleinad · 2023-09-18T19:22:34Z

src/interpolation/details/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp

+    i++;
+    for (; i < q; i++)


Suggested change

i++;

for (; i < q; i++)

for (int i = p + 1; i < q; i++)

masterleinad · 2023-09-18T19:23:13Z

src/interpolation/details/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp

+    i++;
+    for (; i < size; i++)


Suggested change

i++;

for (; i < size; i++)

for (int i = q + 1; i < size; i++)

masterleinad · 2023-09-18T19:29:11Z

src/interpolation/details/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp

+    for (int i = 0; i < size; i++)
+    {
+      auto const u_ip = U(i, p);
+      auto const u_iq = U(i, q);
+      U(i, p) = cos_theta * u_ip + sin_theta * u_iq;
+      U(i, q) = -sin_theta * u_ip + cos_theta * u_iq;
+    }


Would it make sense to use hierarchical parallelism for these for-loops or do we assume size to be small enough that it doesn't matter? Or do we have sufficient parallelism already in the outer loop?

It would make a lot of sense to use hierarchical parallelism here (and in other loops) as well as using scratch pads for the auxiliary matrices. However, because this is part of the MLS algorithm, it might not be a bottleneck and it is best to avoid early optimization.

However, the following should be a logical distribution:

Each team treats a single matrix

Threads/Vectors handle the loops and reductions

Use the team scratch pad for ES and U

Fair enough.

Rombur · 2023-09-19T12:47:09Z

src/interpolation/details/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp

+    // +---+---+              +---+---+
+
+    value_t cos_theta, sin_theta, x, y;
+    if (a == c)


Don't you need to have a tolerance here? What happens if a = 1 and c = 1+1e-15?

A tolerance is not really needed as the algorithm behaves correctly even if the values are very close together. This case is here to avoid a "true" division by zero.

Rombur · 2023-09-19T12:57:11Z

src/interpolation/details/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp

+          symmetricPseudoInverseSVDSerialKernel(A, ES, U);
+        });
+  }
+  else if constexpr (InOutMatrices::rank == 2)


Is there a case where we care about using SVD for a single matrix?

No, at least for now. This was added because the interface of the kernel itself takes only a single matrix. However, I will probably remove it as this might be the source of why cuda+nvcc 11.0.3 fails.

Rombur · 2023-09-19T13:17:41Z

test/ArborX_EnableViewComparison.hpp

+  }
+
+  if (!same_dim_size)
+    return;


I think that you don't need this. BOOST_TEST_REQUIRE does what you want (to check).

Not really, as this will abort the current test and not check the other matrices of the same test.

Rombur · 2023-09-19T13:25:56Z

test/tstInterpDetailsSVD.cpp

+    }
+  makeCase<ExecutionSpace, view_t, double[128][128]>(space, 0, mat, inv, 128);
+
+  // Case for invertible 128x128 matrix


Do you plan to add a test in the future?

If I do find an way to easily build a non-trivial 128x128 matrix and its inverse. Since I added this comment I did not found one so I will remove that command. This test could be extended if I ever find said matrix.

masterleinad · 2023-09-19T21:19:58Z

Windows CI is reporting

LINK : fatal error LNK1104: cannot open file 'D:\a\ArborX\ArborX\build\test\headers_self_contained\Debug\ArborX_HeaderSelfContained_interpolation_details_ArborX_InterpDetailsSymmetricPseudoInverseSVD_hpp.exe' [D:\a\ArborX\ArborX\build\test\headers_self_contained\ArborX_HeaderSelfContained_interpolation_details_ArborX_InterpDetailsSymmetricPseudoInverseSVD_hpp.vcxproj]

and all GPU builds seem to be stuck in ArborX_Test_InterpDetailsSVD.

masterleinad · 2023-09-19T21:23:59Z

src/interpolation/details/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp

+  for (int i = 0; i < size; i++)
+    max_eigen = Kokkos::max(Kokkos::abs(ES(i, i)), max_eigen);
+
+  // We inverse the diagonal of ES, except if "0" is found


Suggested change

// We inverse the diagonal of ES, except if "0" is found

// We invert the diagonal of ES, except if "0" is found

masterleinad · 2023-09-19T21:30:57Z

src/interpolation/details/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp

+                               typename ESMatrix::value_type> &&
+                    std::is_same_v<typename ESMatrix::value_type,
+                                   typename UMatrix::value_type>,
+                "Each input matrix must have the same value type");


Suggested change

"Each input matrix must have the same value type");

"All input matrices must have the same value type");

masterleinad · 2023-09-19T21:34:13Z

src/interpolation/details/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp

+// Pseudo-inverse of symmetric matrices using SVD
+// We must find U, E (diagonal and positive) and V such that A = U.E.V^T
+// We also suppose, as the input, that A is symmetric, so U = SV where S is
+// a sign matrix (only 1 or -1 in the diagonal, 0 elsewhere).


Suggested change

// a sign matrix (only 1 or -1 in the diagonal, 0 elsewhere).

// a sign matrix (only 1 or -1 on the diagonal, 0 elsewhere).

masterleinad · 2023-09-19T21:35:25Z

src/interpolation/details/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp

+    }
+
+  static constexpr value_t epsilon = Kokkos::Experimental::epsilon_v<float>;
+  while (true)


I guess the CI is stuck in this loop?

masterleinad · 2023-09-19T21:36:05Z

src/interpolation/details/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp

+    // | ES(q, p) | ES(q, q) |   | b | c |
+    // +----------+----------+   +---+---+
+
+    // Lets compute x, y and theta such that


Suggested change

// Lets compute x, y and theta such that

// Let's compute x, y and theta such that

masterleinad · 2023-09-19T21:36:25Z

src/interpolation/details/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp

+      y = a + c - x;
+    }
+
+    // Now lets compute the following new values for U and ES


Suggested change

// Now lets compute the following new values for U and ES

// Now let's compute the following new values for U and ES

masterleinad · 2023-09-19T21:36:43Z

src/interpolation/details/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp

+    for (int i = 0; i < size; i++)
+    {
+      auto const u_ip = U(i, p);
+      auto const u_iq = U(i, q);
+      U(i, p) = cos_theta * u_ip + sin_theta * u_iq;
+      U(i, q) = -sin_theta * u_ip + cos_theta * u_iq;
+    }


Fair enough.

masterleinad · 2023-09-19T21:37:06Z

src/interpolation/details/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp

+    }
+  }
+
+  // We compute the max to get a range of the invertible eigen values


Suggested change

// We compute the max to get a range of the invertible eigen values

// We compute the max to get a range of the invertible eigenvalues

mrlag31 · 2023-09-20T13:03:05Z

Windows CI is reporting

LINK : fatal error LNK1104: cannot open file 'D:\a\ArborX\ArborX\build\test\headers_self_contained\Debug\ArborX_HeaderSelfContained_interpolation_details_ArborX_InterpDetailsSymmetricPseudoInverseSVD_hpp.exe' [D:\a\ArborX\ArborX\build\test\headers_self_contained\ArborX_HeaderSelfContained_interpolation_details_ArborX_InterpDetailsSymmetricPseudoInverseSVD_hpp.vcxproj]

I don't know why Windows cannot build this. The logs show no error. I am wondering if the headers must be referenced somewhere else.

and all GPU builds seem to be stuck in ArborX_Test_InterpDetailsSVD.

It seems more likely that, in job 12, only cuda 11.0.3 + clang stopped on that test due to the job time limit (3h) and the others copied the logs of that test.

aprokop · 2023-09-20T17:18:03Z

src/interpolation/details/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp

+        auto const val = Kokkos::abs(mat(i, j) - mat(j, i));
+        auto const ref = Kokkos::abs(mat(i, j));
+        static constexpr value_t epsilon =
+            Kokkos::Experimental::epsilon_v<float>;
+        if (ref == value_t(0) && val > epsilon)
+          return false;
+        if (ref != value_t(0) && val / ref > epsilon)
+          return false;


I'm not sure why the check is using tolerances. I think the input matrix should always be symmetric without tolerances, if constructed correctly. So I would just do

Suggested change

auto const val = Kokkos::abs(mat(i, j) - mat(j, i));

auto const ref = Kokkos::abs(mat(i, j));

static constexpr value_t epsilon =

Kokkos::Experimental::epsilon_v<float>;

if (ref == value_t(0) && val > epsilon)

return false;

if (ref != value_t(0) && val / ref > epsilon)

return false;

if (mat(i, j) != mat(j, i))

return false;

If that's not the case, I would rename the function to include Tol in the name, and provide an epsilon function argument.

I will change to the tolerance-free version, as it forces users to have their matrices properly symmetric.

aprokop · 2023-09-20T17:33:13Z

test/tstInterpDetailsSVD.cpp

+#include "BoostTest_CUDA_clang_workarounds.hpp"
+#include <boost/test/unit_test.hpp>
+
+namespace axid = ArborX::Interpolation::Details;


What does axid stand for?

It stands for ArborX Interpolation Details. It is simply a namespace shortcut. But because it is not used in a lot of places, I will use the full namespace name.

aprokop · 2023-09-20T17:34:19Z

test/tstInterpDetailsSVD.cpp

+  host_view src("src", m, n, n);
+  host_view ref("ref", m, n, n);
+  U inv("inv", m, n, n);


Suggested change

host_view src("src", m, n, n);

host_view ref("ref", m, n, n);

U inv("inv", m, n, n);

host_view src("Testing::src", m, n, n);

host_view ref("Testing::ref", m, n, n);

U inv("Testing::inv", m, n, n);

aprokop · 2023-09-20T17:41:40Z

It seems more likely that, in job 12, only cuda 11.0.3 + clang stopped on that test due to the job time limit (3h) and the others copied the logs of that test.

I'm confused about what's going on with the testing.

mrlag31 · 2023-09-22T16:03:11Z

The windows CI fails because of the location of the added header (two folder deep from src). Jenkins passes with no issues.

masterleinad · 2023-09-22T17:01:25Z

The windows CI fails because of the location of the added header (two folder deep from src). Jenkins passes with no issues.

There are still some warnings to fix, see https://cloud.cees.ornl.gov/jenkins-ci/job/ArborX/job/PR-950/19/gcc/.

masterleinad · 2023-09-22T20:20:55Z

We need to deal with the failing Windows CI one way or the other. As I said earlier, I think it's OK to disable the header self-containment test for that CI build and open an issue for it.

mrlag31 · 2023-09-25T13:03:52Z

I have removed the self-containment test from the Windows CI and it does pass, and made an issue #953. CUDA 11.5 + MPI pipeline did not start but all other succeeded.

masterleinad

Looks good enough.

mrlag31 added 4 commits September 11, 2023 09:08

SVD from arborx#946

a9b9eea

const and better loops

6f67dd3

MD view comparison

8acf635

Execution space on allocations, layout agnostics test

98c2dd8

aprokop added the enhancement New feature or request label Sep 11, 2023

dalg24 reviewed Sep 12, 2023

View reviewed changes

mrlag31 added 3 commits September 12, 2023 11:28

Extra assertions and better readability

a4f7326

More verbose names and extra empty test

79580e3

Better multi dimensional comparison

797fd65

aprokop reviewed Sep 13, 2023

View reviewed changes

mrlag31 added 2 commits September 13, 2023 15:26

Renaming and switching size_t to int

d7c4ee5

Extra symmetric check

d3b6e8d

dalg24 reviewed Sep 14, 2023

View reviewed changes

mrlag31 marked this pull request as ready for review September 15, 2023 13:50

boost test compliance and renaming pre-conditions

3d105a1

mrlag31 force-pushed the symmetric-pinv-svd branch from 0618fcf to 3d105a1 Compare September 15, 2023 13:59

mrlag31 added 2 commits September 15, 2023 11:07

symmetric evaluation in lambda

20f906a

small fixes and slighly better readability

0c5bb0d

mrlag31 force-pushed the symmetric-pinv-svd branch from b8c9556 to 0c5bb0d Compare September 15, 2023 19:18

masterleinad reviewed Sep 18, 2023

View reviewed changes

Better readability for loops

e2e5838

Rombur reviewed Sep 19, 2023

View reviewed changes

Removing handling of single matrix

4568d88

masterleinad reviewed Sep 19, 2023

View reviewed changes

Typos and comments fix

8af21a9

aprokop reviewed Sep 20, 2023

View reviewed changes

Testing labels and stricter symmetric test

61fff17

mrlag31 force-pushed the symmetric-pinv-svd branch from 6bde40a to 61fff17 Compare September 20, 2023 19:19

Compilation and clang-tidy warning handling

8098fc0

mrlag31 force-pushed the symmetric-pinv-svd branch from 25f374f to 8098fc0 Compare September 22, 2023 14:06

mrlag31 added 2 commits September 22, 2023 14:07

Extra warning suppression

2dfac38

Better "makeCase" template and type edits

e68f83b

deactivating header self-containment for windows ci

4edd575

mrlag31 mentioned this pull request Sep 25, 2023

Windows CI fails on header self-containment if the generated program name is too long #953

Open

masterleinad approved these changes Sep 25, 2023

View reviewed changes

aprokop approved these changes Sep 27, 2023

View reviewed changes

aprokop merged commit ae00f3d into arborx:master Sep 27, 2023
1 check passed

mrlag31 mentioned this pull request Sep 27, 2023

Compact radial basis functions and generic polynomial basis / Utility for moving least squares #954

Merged

		value_t const u = (2 * b) / (a - c);
		value_t const v = 1 / Kokkos::sqrt(u * u + 1);

	// We inverse the diagonal of ES, except if "0" is found
	// We invert the diagonal of ES, except if "0" is found

	"Each input matrix must have the same value type");
	"All input matrices must have the same value type");

	// a sign matrix (only 1 or -1 in the diagonal, 0 elsewhere).
	// a sign matrix (only 1 or -1 on the diagonal, 0 elsewhere).

	// Lets compute x, y and theta such that
	// Let's compute x, y and theta such that

	// Now lets compute the following new values for U and ES
	// Now let's compute the following new values for U and ES

	// We compute the max to get a range of the invertible eigen values
	// We compute the max to get a range of the invertible eigenvalues

Pseudo-inverse of symmetric matrices using SVD / Utility for moving least squares #950

Pseudo-inverse of symmetric matrices using SVD / Utility for moving least squares #950

Conversation

mrlag31 commented Sep 11, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mrlag31 Sep 19, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

masterleinad commented Sep 19, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mrlag31 commented Sep 20, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mrlag31 Sep 20, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

aprokop commented Sep 20, 2023

mrlag31 commented Sep 22, 2023

masterleinad commented Sep 22, 2023

masterleinad commented Sep 22, 2023

mrlag31 commented Sep 25, 2023 • edited Loading

masterleinad left a comment

Choose a reason for hiding this comment

mrlag31 commented Sep 11, 2023 •

edited

Loading

mrlag31 Sep 19, 2023 •

edited

Loading

mrlag31 commented Sep 20, 2023 •

edited

Loading

mrlag31 Sep 20, 2023 •

edited

Loading

mrlag31 commented Sep 25, 2023 •

edited

Loading