Implement af::pinverse() #2279

mark-poscablo · 2018-08-14T19:14:27Z

SVD-based Moore-Penrose pseudoinverse.

This has the same type restrictions as af::inverse(), but lifts the dims restriction of requiring the array to be square. The third dimension is also allowed for batching, although this uses naive batching for now (loop through and successively process each slice along the third dimension). No backend-specific code were added, since all the necessary pieces for the SVD-based approach are already implemented in arrayfire (thus all the code is in src/api/c/pinverse.cpp).

The tests include checking if all the Moore-Penrose conditions hold (see https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_inverse#Definition), for all the the appropriate types (f32, f64, c32, c64).

The CI tests might fail initially until the corresponding PR on arrayfire-data is approved and I update the git submodule reference here.

This addresses #2074, #2077, and #2143 (this last one indirectly).

mark-poscablo · 2018-08-14T19:58:28Z

src/api/c/pinverse.cpp

+    // sVec produced by svd() has minimal dim length (no extra zeroes).
+    // Thus s+ produced by diagCreate() will have minimal dims as well,
+    // and v could have an extra dim0 or u* could have an extra dim1
+    if (v.dims()[1] > sPinvCast.dims()[0]) {


By the way, I'm wondering if this is the fastest implementation of what I'm trying to do in this section. At first, actually, I thought of avoiding creating new arrays (createSubArray()) since this might be expensive especially when the arrays are big. Instead, a reference to only the needed section of the arrays might be better to use - although this could also have the disadvantage of making the data access pattern jump (and thus induce cache misses) near the cropped-out section of the arrays during matmul().

Another way could be padding sVec with 0's instead to match v and u*'s dims for matmul. This creates a new array once (and a smaller one compared to the other arrays), but will waste multiplication and adding with 0's during matmul.

mark-poscablo · 2018-08-14T21:44:36Z

I've been benchmarking this vs just af::svd() and af::svdInPlace(), and I'm thinking of adding an in-place option for pinverse() too, especially since af::svdInPlace() is actually really fast on CUDA (this could be on another PR though). However, it has a couple of issues - I'll write about it in a separate issue.

syurkevi · 2018-08-15T20:42:11Z

src/api/c/pinverse.cpp

+    Array<T> v = transpose(vT, true);
+
+    // Round down small values to zero to avoid large reciprocals later
+    Array<Tr> eps = createValueArray<Tr>(sVec.dims(), scalar<Tr>(1e-6));


Change tolerance to match t = ε⋅max(m,n)⋅max(Σ)

Implemented that relative tolerance calculation but epsilon will be 1e-6 by default instead of machine epsilon. There's an issue with SVD that prevents me from lowering this for float. I should probably write about that issue, but basically if I SVD a given array, make one of the singular values really small (1e-12 perhaps), multiply them back together, and SVD it again, the really small singular value isn't so small anymore (becomes something on the order of 1e-5). Thus I think our SVD cannot produce really small singular values and so I set the default tolerance that high.

syurkevi · 2018-08-15T22:01:34Z

src/api/c/pinverse.cpp

+    }
+    if (uT.dims()[0] > sPinvCast.dims()[1]) {
+        std::vector<af_seq> seqs = {
+            {0., static_cast<double>(sPinvCast.dims()[1]), 1.},


Inclusive! Make exclusive.

syurkevi · 2018-08-21T15:49:46Z

include/af/lapack.h

@@ -200,6 +200,25 @@ namespace af
    */
    AFAPI array inverse(const array &in, const matProp options = AF_MAT_NONE);

+    /**


#if AF_API_VERSION >= 37

Adding this now

syurkevi · 2018-08-21T15:53:52Z

src/api/c/pinverse.cpp

+#include <svd.hpp>
+#include <transpose.hpp>
+
+using af::dim4;


using std::vector
using std::swap

Adding this now too

syurkevi · 2018-08-21T16:10:58Z

src/api/c/pinverse.cpp

+        }
+
+        double validTol = tol;
+        if (validTol < 0.) {


ARG_ASSERT?
User should probably know their tolerance is bad...

pavanky · 2018-08-25T06:26:26Z

@mark-poscablo @syurkevi Try to implement this for non-square matrices inside inverse rather than adding a new function.

umar456 · 2018-08-25T17:35:33Z

@pavanky The extra parameters are not necessary for the regular inverse. We can make this part of the regular inverse and pick good defaults so that inverse will also perform a pseudo inverse if the matrices are not square.

mark-poscablo · 2018-08-27T01:51:19Z

@pavanky I was actually planning to do this at first, but @syurkevi advised me early on that the user most probably wants to explicitly differentiate the pseudoinverse from the regular inverse. But true, as @umar456 said, we could add it as a fallback for regular inverse in case the matrices aren't square - if we do this though I think we should add a boolean parameter to af::inverse that decides whether to allow fallback to pseudoinverse or not, setting it default to false to keep the original API's behavior. Reason is that the user might still want to be alerted that the input matrix isn't square, versus quietly selecting pseudoinverse if that's the case. At least with that extra parameter, the user will have to intentionally set it to true if they really want that quiet selection of pseudoinverse.

umar456

This is looking great. I have made a couple of comments regarding the documentation and style.

umar456 · 2018-08-27T16:54:52Z

docs/details/lapack.dox

+
+\brief Pseudo-invert a matrix
+
+This function calculates the Moore-Penrose pseudoinverse of a matrix \f$A\f$, using \ref af::svd at its core. If \f$A\f$ is of size **M x N**, then its pseudoinverse \f$A^+\f$ will be of size **N x M**.


Please limit lines to 80 characters per line.

Yup I followed suit on the other functions' existing docs, thinking that maybe doxygen somehow requires everything on a single line or else an unnecessary line break will occur. I'll change it to 80 chars per line though

Alright it's good, it automatically adds a space for every new line i add on the .dox file.

umar456 · 2018-08-27T17:03:12Z

include/af/lapack.h

+       \param[in] options determining various properties of matrix \p in
+       \returns \p x, the inverse of the input matrix
+
+       \note \p tol is not the actual lower threshold, but it is passed in as a parameter to the calculation of the actual threshold relative to the shape and contents of \p in.


80 characters per line. You can tab to the first character of the comment to make it more readable:

\note \p tol is not the actual lower threshold, but it is passed in as a parameter to the calculation of the actual threshold relative to the shape and contents of \p in.

Same case as above for this file, just added line breaks now (to pinverse and af_pinverse)

umar456 · 2018-08-27T17:03:42Z

src/api/c/pinverse.cpp

@@ -0,0 +1,170 @@
+/*******************************************************
+ * Copyright (c) 2014, ArrayFire


This is a new file so you should set this to 2018

Makes sense, I just changed this one and test/pinverse.cpp

umar456 · 2018-08-27T17:04:14Z

src/api/c/pinverse.cpp

+Array<T> pinverseSvd(const Array<T> &in, const double tol)
+{
+    // Moore-Penrose Pseudoinverse
+


Nit: Extra line

Here I intentionally put an extra line because the comment doesn't apply only to the code block right below it, but rather the whole function. Maybe I should put this comment above the signature instead eh?

It would be better to apply that comment to the whole function. Maybe move it above the function

umar456 · 2018-08-27T17:07:56Z

src/api/c/pinverse.cpp

+        }
+
+        ARG_ASSERT(1, i_info.isFloating()); // Only floating and complex types
+        ARG_ASSERT(2, tol >= 0.); // Only floating and complex types


The comment is incorrect here.

Whoops copy paste. Just changed the comment.

umar456 · 2018-08-27T17:20:40Z

test/pinverse.cpp

+
+// Test Moore-Penrose conditions
+// See https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_inverse#Definition
+


Nit: Extra line

Same reasoning as in api/c/pinverse.cpp but I just removed the extra line here and clarified that this comment just applies to the first 4 following tests

umar456 · 2018-08-27T17:25:03Z

docs/details/lapack.dox

+
+This function calculates the Moore-Penrose pseudoinverse of a matrix \f$A\f$, using \ref af::svd at its core. If \f$A\f$ is of size **M x N**, then its pseudoinverse \f$A^+\f$ will be of size **N x M**.
+
+This calculation can be batched if the input array is three-dimensional (**M x N x P**). Each **M x N** slice along the third dimension will have its own pseudoinverse, for a total of **P** pseudoinverses in the output array (**N x M x P**).


Instead of making the M x N x P statement bold perhaps use \f$M \times N \times P$\f.

$M \times N \times P$

Yup I was on the fence about this one, since other functions' docs (like af::svd) just usually put bold. But it's probably better to be consistent within a doc page. I'll change this.

umar456 · 2018-08-27T17:27:25Z

include/af/lapack.h

+       \param[in] in is the input matrix
+       \param[in] tol defines the lower threshold for singular values from SVD
+       \param[in] options determining various properties of matrix \p in
+       \returns \p x, the inverse of the input matrix


I don't see an x value in the function signature. just use
\returns the inverse of the input matrix

True, let me take it out.

umar456 · 2018-08-27T17:27:45Z

include/af/lapack.h

+
+       \param[in] in is the input matrix
+       \param[in] tol defines the lower threshold for singular values from SVD
+       \param[in] options determining various properties of matrix \p in


Singular option.

Perhaps state that this must be AF_MAT_NONE

\param[in] option must be AF_MAT_NONE. For future use.

Yup just changed it.

src/api/c/pinverse.cpp

umar456 · 2018-08-27T21:03:33Z

src/api/c/pinverse.cpp

+                {0., static_cast<double>(inArray.dims()[1] - 1), 1.},
+                {static_cast<double>(i), static_cast<double>(i), 1.}
+            };
+            Array<T> inSlice = createSubArray<T>(inArray, seqs);


Can you make sure that this function is indeed making a subarray and not allocating and performing a copy? If it's copying an array then we should move this operation into the svd function and then iterate over the arrays there. You can just create a TODO if we are doing a copy and we can fix it at a later time. I suspect it is working as expected but I just want to make sure

…rrectly write into subarrays (arrayfire#2279)

9prady9 · 2018-09-19T04:05:13Z

src/api/c/pinverse.cpp

+    Array<T> vT = createValueArray<T>(dim4(N, N, P, Q), scalar<T>(0));
+    Array<Tr> sVec = createValueArray<Tr>(dim4(min(M, N), 1, P, Q), scalar<Tr>(0));
+    for (uint j = 0; j < Q; ++j) {
+        for (uint i = 0; i < P; ++i) {


Not sure if it will effect performance, but this double-loop can be merged into a single loop in this case.

src/api/c/pinverse.cpp

…rrectly write into subarrays (arrayfire#2279) (cherry picked from commit 6fc326f)

…rrectly write into subarrays (#2279) (cherry picked from commit 6fc326f)

mark-poscablo requested a review from umar456 August 14, 2018 19:15

mark-poscablo commented Aug 14, 2018

View reviewed changes

mark-poscablo mentioned this pull request Aug 14, 2018

svdInPlace() does not accept arrays with dim0 > dim1 #2282

Closed

mark-poscablo force-pushed the pinverse branch 3 times, most recently from b75a6f4 to 89c357d Compare August 17, 2018 23:59

mark-poscablo added the feature label Aug 24, 2018

mark-poscablo modified the milestone: v3.7.0 Aug 24, 2018

syurkevi reviewed Aug 24, 2018

View reviewed changes

umar456 requested changes Aug 27, 2018

View reviewed changes

umar456 reviewed Aug 27, 2018

View reviewed changes

mark-poscablo force-pushed the pinverse branch from 2964698 to 8d1a81c Compare August 30, 2018 18:55

mark-poscablo force-pushed the pinverse branch 2 times, most recently from 1ed87a2 to d592c34 Compare September 10, 2018 21:26

mark-poscablo added 2 commits September 12, 2018 15:26

Added pinverse (arrayfire#2279)

70a6a89

svd OpenCL: Use buffer map/unmap instead of read/write in order to co…

cc90985

…rrectly write into subarrays (arrayfire#2279)

mark-poscablo force-pushed the pinverse branch from d592c34 to cc90985 Compare September 12, 2018 19:40

umar456 approved these changes Sep 19, 2018

View reviewed changes

9prady9 reviewed Sep 19, 2018

View reviewed changes

src/api/c/pinverse.cpp Show resolved Hide resolved

umar456 merged commit 6fc326f into arrayfire:master Sep 20, 2018

umar456 pushed a commit that referenced this pull request Sep 20, 2018

Added pinverse (#2279)

b4f9230

umar456 pushed a commit to umar456/arrayfire that referenced this pull request Nov 2, 2018

svd OpenCL: Use buffer map/unmap instead of read/write in order to co…

612bc1d

…rrectly write into subarrays (arrayfire#2279) (cherry picked from commit 6fc326f)

umar456 pushed a commit to umar456/arrayfire that referenced this pull request Nov 2, 2018

svd OpenCL: Use buffer map/unmap instead of read/write in order to co…

795d291

…rrectly write into subarrays (arrayfire#2279) (cherry picked from commit 6fc326f)

umar456 pushed a commit to umar456/arrayfire that referenced this pull request Nov 3, 2018

svd OpenCL: Use buffer map/unmap instead of read/write in order to co…

8447030

…rrectly write into subarrays (arrayfire#2279) (cherry picked from commit 6fc326f)

9prady9 pushed a commit that referenced this pull request Nov 3, 2018

svd OpenCL: Use buffer map/unmap instead of read/write in order to co…

af799fc

…rrectly write into subarrays (#2279) (cherry picked from commit 6fc326f)


		\brief Pseudo-invert a matrix

		This function calculates the Moore-Penrose pseudoinverse of a matrix \f$A\f$, using \ref af::svd at its core. If \f$A\f$ is of size M x N, then its pseudoinverse \f$A^+\f$ will be of size N x M.

		@@ -0,0 +1,170 @@
		/*******************************************************
		* Copyright (c) 2014, ArrayFire


		// Test Moore-Penrose conditions
		// See https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_inverse#Definition


		This function calculates the Moore-Penrose pseudoinverse of a matrix \f$A\f$, using \ref af::svd at its core. If \f$A\f$ is of size M x N, then its pseudoinverse \f$A^+\f$ will be of size N x M.

		This calculation can be batched if the input array is three-dimensional (M x N x P). Each M x N slice along the third dimension will have its own pseudoinverse, for a total of P pseudoinverses in the output array (N x M x P).

Implement af::pinverse() #2279

Implement af::pinverse() #2279

Uh oh!

Conversation

mark-poscablo commented Aug 14, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

mark-poscablo Aug 14, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mark-poscablo commented Aug 14, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

pavanky commented Aug 25, 2018

Uh oh!

umar456 commented Aug 25, 2018

Uh oh!

mark-poscablo commented Aug 27, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

umar456 left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mark-poscablo commented Aug 14, 2018 •

edited

Loading

mark-poscablo Aug 14, 2018 •

edited

Loading

mark-poscablo commented Aug 14, 2018 •

edited

Loading

mark-poscablo commented Aug 27, 2018 •

edited

Loading