Replace lapack with linalg in MahalanobisDistance #4085

vinx13 · 2018-01-16T13:20:07Z

This PR replaces lapack usage with linalg in MahalanobisDistance.
A unit test of MahalanobisDistance is added (reference computed with Matlab)

karlnapf · 2018-01-16T13:21:14Z

Great!

karlnapf

Great, thanks for that.

Is the mahalanobis covered by an integration test? would be something that matches integration_meta_cpp*mahalanobis*

I think we could also use linalg in ridge regression, rather than pure eigen3

karlnapf · 2018-01-16T13:21:38Z

cmake/FindMetaExamples.cmake

@@ -19,7 +19,6 @@ function(get_excluded_meta_examples)
    read_config_header(${shogun_INCLUDE_DIR}/shogun/lib/config.h)
    is_set(${SHOGUN_CONFIG} HAVE_NLOPT)
    is_set(${SHOGUN_CONFIG} USE_GPL_SHOGUN)
-    is_set(${SHOGUN_CONFIG} HAVE_LAPACK)


we want to keep that, there might be other examples

karlnapf · 2018-01-16T13:22:19Z

cmake/FindMetaExamples.cmake

@@ -37,13 +36,6 @@ function(get_excluded_meta_examples)
 			)
 	ENDIF()

-	IF(NOT HAVE_LAPACK)
-		LIST(APPEND EXCLUDED_META_EXAMPLES


please keep the block, so we can add things there that pop up in the future, just don't append anything

karlnapf · 2018-01-16T13:25:24Z

src/shogun/regression/LinearRidgeRegression.cpp

@@ -8,7 +8,6 @@
 */
 #include <shogun/lib/config.h>

-#ifdef HAVE_LAPACK


ah so this guard was just outdated?

I think we should use linalg rather than eigen3 directly inside though

Yes, just outdated. I will look into that to find the possibility of using pure linalg.

I found eigen_kernel_matrix.diagonal() += eigen_tau;. We are not able to implement adding a vector to a diagonal with SGMatrix or linalg. Any suggestions?

leave it as it is for now (you can already change the rest)
Then you could create this method in linalg, and then send another PR

Or, you could just add a loop for now.

vinx13 · 2018-01-16T13:44:16Z

both MahalanobisDistance and LinearRidgeRegression are covered by integration tests.

karlnapf · 2018-01-17T17:55:52Z

I think the LAPACK guard was still necessary.
BTW if you remove such guards, make sure that Shogun compiles if the library that you guarded is not available. Really important.

The windows build doesnt have lapack installed, and here you go:

[00:29:53]   C:\projects\shogun\src\shogun\distance\MahalanobisDistance.cpp(56): error C2039: 'inverse': is not a member of 'shogun::SGMatrix<float64_t>' [C:\projects\shogun\build\src\shogun\libshogun.vcxproj]
[00:29:53]   C:\projects\shogun\src\shogun\distance\MahalanobisDistance.cpp(56): error C3861: 'inverse': identifier not found [C:\projects\shogun\build\src\shogun

karlnapf · 2018-01-17T17:59:50Z

It is this line:
https://github.com/shogun-toolbox/shogun/blob/develop/src/shogun/distance/MahalanobisDistance.cpp#L58

But actually this code is bad generally (not your fault). One should never ever store matrix inverse. Instead, it should store a cholesky factorization of the covariance matrix. Then, rather than multiplying the inverse covariance matrix to a vector in the CMahalanobisDistance::compute method, you should do a cholesky (triangular) solve. See (https://github.com/shogun-toolbox/shogun/blob/develop/src/shogun/classifier/LDA.cpp#L180 for an example how to do this in linalg

vinx13 · 2018-01-18T06:39:11Z

@karlnapf Since cov matrix is PSD, we should use LDLT (linalg is using LLT) for Cholesky decomposition, right? Shall I add it to linalg?

karlnapf · 2018-01-18T10:03:45Z

The matrix is PSD indeed, so LDLT would be more stable.
https://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html

I think adding an argument to cholesky would be good, so that the factorization can operate in both modes, something like or use_ldlt which by default is false. Also make sure to update the docs appropriately so it is clear what is used.

Good catch!

vinx13 · 2018-01-18T15:03:40Z

Seems simply adding a use_ldlt option will not work. In LDLT, we want both L and D. One solution is to overload cholesky_factor as SGMatrix<T> cholesky_factor(const SGMatrix<T>& A, SGMatrix<T> &D). Since D is diagonal, alternatively we could use a SGVector. Besides, we may need to add a LDLT solver.

karlnapf · 2018-01-18T18:34:09Z

You are absolutely right. In fact I now remember a discussion we previously had about this some years back. I think a clean solution would be to expose a new method and a few solver method as well. Everything else is too messy.

OXPHOS · 2018-01-20T15:10:38Z

We did have some discussion years ago. The best way seems to be having separate methods for LLT and LDLT. And since now the linalg methods are split into several files, having more than one cholesky decomposition or solver would be less messy...

vinx13 · 2018-01-20T15:25:05Z

I see. I'll try to implement that.

karlnapf · 2018-02-05T09:46:04Z

src/shogun/distance/MahalanobisDistance.h

+	 * @see <a href="http://en.wikipedia.org/wiki/Mahalanobis_distance">
+	 * Wikipedia: Mahalanobis Distance</a>
+	 */
+	class CMahalanobisDistance : public CRealDistance


BTW why this reformatting?
Pls don't put this into a PR that changes small bits as this one. Makes everything much harder to review

Seems the style checker changed other codes than mine. I'm using command like git clang-format-3.8 --commit xxx , is it right ?

karlnapf · 2018-02-05T09:46:28Z

src/shogun/regression/LinearRidgeRegression.cpp

-    eigen_kernel_matrix = eigen_feats_matrix*eigen_feats_matrix.transpose();
-
-    eigen_kernel_matrix.diagonal() += eigen_tau;
+	linalg::matrix_prod(feats_matrix, feats_matrix, kernel_matrix, false, true);


nice much cleaner

karlnapf · 2018-02-05T09:46:45Z

src/shogun/regression/LinearRidgeRegression.cpp

-    	return false;
-    }
-    eigen_y = llt.solve(eigen_y);
+	auto decomposition = linalg::cholesky_factor(kernel_matrix);


why did you change the indentation? Is this the style cheker?

Just noticed that the original code is indented with spaces.

karlnapf · 2018-02-05T09:47:37Z

tests/unit/distance/MahalanobisDistance_unittest.cc

@@ -0,0 +1,58 @@
+/*
+ * Copyright (c) The Shogun Machine Learning Toolbox
+ * Written (w) 2018 Wuwei Lin


We changed the license header recently, to something much shorter. Check other files for the template ...

karlnapf · 2018-02-05T09:49:12Z

Could you pls split the PR into two: one for regression and one for the distance.
You still haven't checked whether shogun compiles with the guard remove and no lapack installed.
Pls don't remove guards without checking whether things still work: https://ci.appveyor.com/project/vigsterkr/shogun/build/1873#L6478

vinx13 · 2018-02-05T10:30:52Z

@karlnapf Created a new PR for LinearRidgeRegression (#4149 ). I will wait until #4130 be merged before I can continue with the refactor of MahalanobisDistance. Can I force update this PR now to drop commits with regression ?

karlnapf · 2018-02-05T13:42:43Z

yep, you can just force update this one now

karlnapf · 2018-02-05T13:43:03Z

But you will need to fix the "inverse" problem (see above) first.

vigsterkr · 2018-02-06T08:24:58Z

src/shogun/distance/MahalanobisDistance.cpp


 	if ( l == r)
 	{
 		mean = ((CDenseFeatures<float64_t>*) l)->get_mean();
-		icov  = ((CDenseFeatures<float64_t>*) l)->get_cov();
+		cov = ((CDenseFeatures<float64_t>*)l)->get_cov();


while you are changing these lines, could you please replace these explicit casts to dynamic_cast at least plz? :)
or even better do a REQUIRED assertation for the type

vigsterkr · 2018-02-06T08:25:07Z

src/shogun/distance/MahalanobisDistance.cpp

 	}
 	else
 	{
 		mean = ((CDenseFeatures<float64_t>*)l)
 		           ->compute_mean((CDotFeatures*)lhs, (CDotFeatures*)rhs);
-		icov = CDotFeatures::compute_cov((CDotFeatures*) lhs, (CDotFeatures*) rhs);
+		cov = CDotFeatures::compute_cov((CDotFeatures*)lhs, (CDotFeatures*)rhs);


another casting

vigsterkr · 2018-02-06T08:25:23Z

src/shogun/distance/MahalanobisDistance.cpp

-
-	float64_t result = cblas_ddot(v.vlen, v.vector, 1, diff.vector, 1);
+	auto v = linalg::ldlt_solver(chol_cov_L, chol_cov_d, chol_cov_p, diff);
+	float64_t result = linalg::dot(v, diff);


this could be auto or? :)

karlnapf · 2018-02-06T10:45:06Z

src/shogun/distance/MahalanobisDistance.cpp

@@ -39,20 +37,43 @@ bool CMahalanobisDistance::init(CFeatures* l, CFeatures* r)
 {
 	CRealDistance::init(l, r);


this call covers all the checks, so we don't need them again in this call (just convolutes the code)
@vigsterkr both REAL and DENSE is checked. That was the whole point of having the base class here iirc

I would do a cast after the call and add a comment

i would actually change the design of this and have instead of init returning a boolean the actual typecasted values available... as i dont see really the point of an assertation failure that will just throw an exception and in case of success have a true :)
this function call should actually tell us that yes the types are cool and here you go this is what you wanted anyways

@vinx13 you can ignore most of this... just make sure that you do an ASSERT(CRealDistance::init(l, r)); and add a //FIXME: above with reference to the above comment

@vigsterkr I mean this whole type system in the features is messed up. I would just rely on dynamic_cast everywhere. The bool is also pointless. And the enums are even worse (too many types, I would make those strings or even better drop them)

vigsterkr · 2018-02-06T10:45:52Z

src/shogun/distance/MahalanobisDistance.cpp

+	    "Right hand side (was %s) has to be CDenseFeatures instance.\n",
+	    rhs->get_name());
+
+	REQUIRE(


and actually this is implicitly tested in CRealDistance::init
https://github.com/vinx13/shogun/blob/39e7dd2f390d13e665678be23b0b0221f1212033/src/shogun/distance/RealDistance.h#L34

haha I wrote the same thing above ;) sorry didnt see

karlnapf · 2018-02-06T10:46:05Z

src/shogun/distance/MahalanobisDistance.cpp

@@ -63,9 +84,10 @@ void CMahalanobisDistance::cleanup()

 float64_t CMahalanobisDistance::compute(int32_t idx_a, int32_t idx_b)
 {
+	auto feat_l = static_cast<CDenseFeatures<float64_t>*>(lhs);


this is reliable due to old code design principles by sören

vigsterkr · 2018-02-06T10:47:16Z

src/shogun/distance/MahalanobisDistance.cpp

@@ -39,20 +37,43 @@ bool CMahalanobisDistance::init(CFeatures* l, CFeatures* r)
 {
 	CRealDistance::init(l, r);


instead of the REQUIRED sections it should be just check the return type of CRealDistance::init(l, r) as it checks for everything you need want....

karlnapf · 2018-02-06T10:47:44Z

src/shogun/distance/MahalanobisDistance.cpp

-		diff.vlen, diff.vector, 1, 0.0, v.vector, 1);
-
-	float64_t result = cblas_ddot(v.vlen, v.vector, 1, diff.vector, 1);
+	auto v = linalg::ldlt_solver(chol_cov_L, chol_cov_d, chol_cov_p, diff);


yeah this is much much from a numerical perspective!

vinx13 · 2018-02-18T17:48:13Z

@karlnapf could you plz take a look

karlnapf · 2018-02-19T17:11:29Z

src/shogun/distance/MahalanobisDistance.cpp

@@ -37,22 +35,31 @@ CMahalanobisDistance::~CMahalanobisDistance()

 bool CMahalanobisDistance::init(CFeatures* l, CFeatures* r)
 {
-	CRealDistance::init(l, r);
+	// FIXME: See comments in


* Fix broken link in doc * Use linalg in MahalanobisDistance * Apply autoformatter * Remove MehalanobisDistance from LAPACK guarded list * Add lapack guard back * Update license header style * Save ldlt instead of inverse of cov * Replace explict casts and add type check * Remove unnecessary type checks

karlnapf requested changes Jan 16, 2018

View reviewed changes

vinx13 mentioned this pull request Jan 18, 2018

Add linalg API to allow adding a vector to diagonal #4087

Merged

vinx13 mentioned this pull request Jan 28, 2018

Add LDLT Cholesky decomposition #4130

Merged

vinx13 force-pushed the feature/use_linalg branch from d2ad7c8 to b2206e6 Compare February 5, 2018 07:17

karlnapf reviewed Feb 5, 2018

View reviewed changes

vinx13 force-pushed the feature/use_linalg branch from a2211df to caa2902 Compare February 5, 2018 13:59

vinx13 changed the title ~~Replace lapack with linalg~~ Replace lapack with linalg in MahalanobisDistance Feb 5, 2018

vinx13 added 6 commits February 6, 2018 14:32

Fix broken link in doc

32a0046

Use linalg in MahalanobisDistance

eadeadc

Apply autoformatter

0b6e504

Remove MehalanobisDistance from LAPACK guarded list

b9382bf

Add lapack guard back

66f133f

Update license header style

fb42c9a

vinx13 force-pushed the feature/use_linalg branch from caa2902 to fb42c9a Compare February 6, 2018 06:51

Save ldlt instead of inverse of cov

e93e846

vigsterkr requested changes Feb 6, 2018

View reviewed changes

Replace explict casts and add type check

39e7dd2

vigsterkr approved these changes Feb 6, 2018

View reviewed changes

shogun-toolbox deleted a comment from vigsterkr Feb 6, 2018

karlnapf reviewed Feb 6, 2018

View reviewed changes

vigsterkr reviewed Feb 6, 2018

View reviewed changes

karlnapf reviewed Feb 6, 2018

View reviewed changes

vigsterkr reviewed Feb 6, 2018

View reviewed changes

karlnapf reviewed Feb 6, 2018

View reviewed changes

Remove unnecessary type checks

4a4dffa

karlnapf reviewed Feb 19, 2018

View reviewed changes

karlnapf approved these changes Feb 19, 2018

View reviewed changes

karlnapf merged commit bb72427 into shogun-toolbox:develop Feb 19, 2018

		@@ -39,20 +37,43 @@ bool CMahalanobisDistance::init(CFeatures* l, CFeatures* r)
		{
		CRealDistance::init(l, r);

Replace lapack with linalg in MahalanobisDistance #4085

Replace lapack with linalg in MahalanobisDistance #4085

Conversation

vinx13 commented Jan 16, 2018 • edited

karlnapf commented Jan 16, 2018

karlnapf left a 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

vinx13 commented Jan 16, 2018

karlnapf commented Jan 17, 2018

karlnapf commented Jan 17, 2018 • edited

vinx13 commented Jan 18, 2018

karlnapf commented Jan 18, 2018

vinx13 commented Jan 18, 2018

karlnapf commented Jan 18, 2018

OXPHOS commented Jan 20, 2018

vinx13 commented Jan 20, 2018

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

karlnapf commented Feb 5, 2018

vinx13 commented Feb 5, 2018 • edited

karlnapf commented Feb 5, 2018

karlnapf commented Feb 5, 2018

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

vinx13 commented Feb 18, 2018

Choose a reason for hiding this comment

vinx13 commented Jan 16, 2018 •

edited

karlnapf commented Jan 17, 2018 •

edited

vinx13 commented Feb 5, 2018 •

edited