Add Cholesky Decompostion

[UmashankarTriforce]

This PR aims to add Cholesky Decomposition as per the issue #1448

Route taken to complete this feature -

• Implement function
• Write unit test cases
• Document the function

[toslunar]

Why don't you use numpy.linalg.cholesky?

[UmashankarTriforce]

How is cupy.linalg.cholesky implemented? I'm using cupy.linalg.cholesky for backward computation. I checked the difference between numpy and scipy implementation of cholesky.

scipy.linalg.cholesky is giving you the upper-triangular decomposition by default, whereas np.linalg.cholesky is giving you the lower-triangular version.

https://stackoverflow.com/questions/16699163/what-is-the-difference-between-cholesky-in-numpy-and-scipy

[toslunar]

How is cupy.linalg.cholesky implemented?

cupy.linalg.cholesky decomposes a matrix into L @ L.H. The code https://github.com/cupy/cupy/blob/c9c7c4eedcca54e6f53b5dcfc205d4ad94a2b595/cupy/linalg/decomposition.py#L49-L54 has cublas.CUBLAS_FILL_MODE_UPPER, which looks inconsistent but it's not because of the difference of row-major vs column-major.

SciPy-compatible GPU implementation is available in cupyx.scipy namespace, if there is.

I'm using cupy.linalg.cholesky for backward computation.

Inputs to FunctionNode.backward are Variables.

• Use Chainer functions in backward, or
• define class CholeskyGrad(FunctionNode) and leave double-backward unimplemented.

[UmashankarTriforce]

Sorry for the late response!
SciPy-compatible GPU implementation is not there for cholesky. Also I plan to use cupy.linalg.cholesky in forward_gpu and for backwards I plan to take the dot product of L and L.T.conj(). Since chainer.dot is not implemented, I plan to use cupy.dot and wrap the array back to chainer.Variable. What do you suggest?

Edit: I also found that the results of cupy.linalg.cholesky is consistent with numpy's implementation. So I would be switching to numpy for forward_cpu

[toslunar]

Variable-level equivalent code would be

a, = self.get_retrained_inputs()
return chainer.functions.matmul(a, a, transb=True),

[UmashankarTriforce]

Sorry for the late response! I have added changes and also added test patterns
There's an error I'm encountering when I try to run the test
AttributeError: 'Cholesky' object has no attribute 'get_retrained_inputs'

How do I solve this?

[UmashankarTriforce]

When I run tests, in some cases the gradients computed backwards is absurdly more than the numeric gradients. For example -

>           raise AssertionError(f.getvalue())
E           chainer.testing.function_link.FunctionTestError: Parameterized test failed.
E           
E           Base test method: TestCholesky_use_chainerx_false__chainerx_device_None__use_cuda_false__cuda_device_None__use_cudnn_never__cudnn_deterministic_false__autotune_false__cudnn_fast_batch_normalization_false__use_ideep_never.test_backward
E           Test parameters:
E             dtype: <class 'numpy.float32'>
E             shape: (2, 2)
E           
E           
E           (caused by)
E           FunctionTestError: backward is not implemented correctly
E           
E           (caused by)
E           AssertionError: check_backward failed (eps=0.001 atol=0.001 rtol=0.001)
E           inputs[0]:
E           [[7670.438  1473.2218]
E            [1473.2218 6982.2075]]
E           grad_outputs[0]:
E           [[-0.70425826  0.97613996]
E            [ 0.5500437  -0.09768751]]
E           directions[0]:
E           [[-0.7262568  -0.0971007 ]
E            [-0.21805311  0.64465135]]
E           gradients (numeric):  0.0015698603367818734
E           gradients (backward): 23360998.461244226
E           
E           x: numeric gradient, y: backward gradient
E           Not equal to tolerance rtol=0.001, atol=0.001
E           
E           Mismatch: 100%
E           Max absolute difference: 23360998.45967437
E           Max relative difference: 1.
E            x: array(0.00157)
E            y: array(23360998.461244)
E           
E           assert_allclose failed: 
E             shape: () ()
E             dtype: float64 float64
E             i: (0,)
E             x[i]: 0.0015698603367818734
E             y[i]: 23360998.461244226
E             relative error[i]: 0.9999999999328
E             absolute error[i]: 23360998.459674366
E             relative tolerance * |y[i]|: 23360.998461244228
E             absolute tolerance: 0.001
E             total tolerance: 23360.999461244228
E           x: 0.00156986
E           y: 23360998.46124423
``
How do you suggest I fix this issue?

[toslunar]

The scale in the CuPy test is large because integer types are tested. To test backprop, scale=(1e-2, 2.0) looks good to me.

[UmashankarTriforce]

Yes I have made the necessary changes, but the mismatch still exists (its of a lower magnitude now, but more than tolerance value)

[toslunar]

My fix is available at UmashankarTriforce#1. Could you merge it?

[UmashankarTriforce]

Yes its done

[toslunar]

• In order to add the link to the function, could you edit docs/source/reference/functions.rst?
• On the backprop, I changed my mind and let's return a symmetric matrix. Please wait and see Document properties of computed gradients in cholesky and eigh #8312.
• F.cholesky could call chainerx.linalg.cholesky by overriding forward_chainerx:

    def forward_chainerx(self, inputs):
        return chainerx.linalg.cholesky(*inputs),

[UmashankarTriforce]

I have added the points 1 and 3 from your review. I will wait for #8312 and change accordingly

[toslunar]

#8312 is merged. Could you align F.cholesky with chainerx.linalg.cholesky?

[toslunar]

Insert a = 0.5 * (a + a.T) to make the perturbation symmetric, like the test code for chainerx.linalg.cholesky:

chainer/tests/chainerx_tests/unit_tests/routines_tests/test_linalg.py

Lines 504 to 514 in c8e3b4f

	def forward_xp(self, inputs, xp):
	a, = inputs
	if (_numpy_does_not_support_0d_input113 and a.size == 0):
	pytest.skip('Older NumPy versions do not work with empty arrays')
	# Input has to be symmetrized for backward test to work
	a = (a + a.T)/2. + 1e-3 * xp.eye(*self.shape)
	L = xp.linalg.cholesky(a)
	return L,

.

It can be also tested that the returned gradient is symmetric. (Not mandatory in this PR)

[toslunar]

Could you also fix TestCholesky.forward? xp means numpy or cupy (or chainerx) in Chainer.

[UmashankarTriforce]

Yes I have added this fix

[UmashankarTriforce]

I'm extremely sorry for the late response 🙁 : I had my exams in college and hence I wasn't able to work on this. Have a look and let me know if any more modifications are necessary 🙂

[toslunar]

Could you also fix TestCholesky.forward? xp means numpy or cupy (or chainerx) in Chainer.

[toslunar]

LGTM

[toslunar]

Jenkins & flexCI, test this please.

[chainer-ci]

Jenkins CI test (for commit a496ded, target branch master) failed with status FAILURE.

[UmashankarTriforce]

The error seems to be from ChainerX