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
)
[1]:
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:
[2]:
np.linalg.cholesky(np.eye(100))
[2]:
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:
[3]:
# 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
[4]:
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:
[5]:
A = np.eye(1000, dtype=np.float64)
A[:, 2] = 1
A = np.dot(A, A.T)
[6]:
L = cholesky_timed(A)
Elapsed Time: 31.277 ms (1.1e+07 GOP/s)
[ ]: