Cholesky decomposition

A Cholesky decomposition is a useful factorization of Hermitian, positive-definite matrices into the product of a lower triangular matrix \(L\) with its conjugate transpose \(L^{*}\).

Numpy has a function numpy.linalg.cholesky built-in for computing Cholesky decompositions. cuPyNumeric also implements this function, and it can be used as an immediate drop-in replacement.

License

Copyright 2024 NVIDIA Corporation

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

     http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

To get started, import cupynumeric as np (just the same way we would import numpy)

import cupynumeric as np  # instead of numpy

At this point we can call np.linalg.cholesky, exactly how we would with Numpy, but will get the result computed by cuPyNumeric’s cholesky function. Let’s quickly try it out with a simple identitity matrix:

np.linalg.cholesky(np.eye(100))

array([[1., 0., 0., ..., 0., 0., 0.],
       [0., 1., 0., ..., 0., 0., 0.],
       [0., 0., 1., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 1., 0., 0.],
       [0., 0., 0., ..., 0., 1., 0.],
       [0., 0., 0., ..., 0., 0., 1.]])

We’d like to get some information about how well cuPyNumeric’s cholesky function performs. In order to obtain accurate timings, we need to use the time function from legate.timing. Let’s define a helper function cholesky_timed that calls the time function for us, and prints out the results as well:

# Because of Legate's deferred execution model, legate.timing should be used instead
# of standard Python datetime utilities. Python datetime.now would return the time
# a task is *scheduled*, not necessarily the time a task finishes executing.
from legate.timing import time

def cholesky_timed(A):
    start = time()
    result = np.linalg.cholesky(A)
    stop = time()

    n = A.shape[0]
    flops = (n**3) / 3 + 2 * n / 3
    total = (stop - start) / 1000.0
    print(f"Elapsed Time: {total} ms ({(flops / total):0.2} GOP/s)")

    return result

Now we can define a matrix \(A\) to decompose. Let’s make some thing a little more interesting than a plain identity matrix:

A = np.eye(1000, dtype=np.float64)
A[:, 2] = 1
A = np.dot(A, A.T)

L = cholesky_timed(A)

Elapsed Time: 31.277 ms (1.1e+07 GOP/s)