Skip to content

ENH: Add linalg.solve_triangular #927

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom
Merged
mrocklin merged 2 commits into from
Jan 22, 2016
Merged

ENH: Add linalg.solve_triangular #927

mrocklin merged 2 commits into from
Jan 22, 2016

Conversation

@sinhrks
Copy link
Member

@sinhrks sinhrks commented Jan 19, 2016

Added solve_triangular, going to use with #885 to implement solve and inv.

import numpy as np
import dask.array as da

np.random.seed(1)
A = np.random.random_integers(1, 10, (12, 12))
b = np.random.random_integers(1, 10, 12)
A = np.tril(A)

dA = da.from_array(A, (3, 3))
db = da.from_array(b, 3)

res = da.linalg.solve_triangular(dA, db, lower=True)
res.compute()
# array([   1.33333333,    0.33333333,   -0.83333333,   -0.76190476,
#          15.95238095,  -25.97142857,   -6.11428571,   43.34444444,
#         -18.42638889,  -17.75753968,  -12.45119048,  139.27301587])

np.allclose(dA.dot(res), b)
# True

mrocklin
mrocklin reviewed Jan 19, 2016
View reviewed changes
dask/array/linalg.py Outdated
Copy link
Member

@mrocklin mrocklin Jan 19, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can we move functions like this outside of the solve_triangular function? Dynamically created functions are expensive to serialize. This is important for distributed computation.

@mrocklin
Copy link
Member

mrocklin commented Jan 19, 2016

How hard would it be to extend this to the case where b is a tall-and-skinny matrix with chunks like ((k, k, k, k), (k,))?

@mrocklin
Copy link
Member

mrocklin commented Jan 19, 2016

Also, in case you're interested in doing even more work, the symmetric version of all of these algorithms is often desired (e.g. lu -> cholesky). In common applications you can usually assume that A is symmetric positive definite.

@sinhrks
Copy link
Member Author

sinhrks commented Jan 19, 2016

@mrocklin OK, could support "tall-and-skinny matrix".

Because solve_triangular itself doesn't rely on lu or choleskey, we can add an option (or separate func) to switch them for intermediate internal process of solve or inv.

@sinhrks
Copy link
Member Author

sinhrks commented Jan 19, 2016

Ah, found that the logic for chunked 2d array is simpler than I thought. Implemented.

@mrocklin
Copy link
Member

mrocklin commented Jan 19, 2016

Ah, found that the logic for chunked 2d array is simpler than I thought. Implemented.

Excellent :)

@mrocklin
Copy link
Member

mrocklin commented Jan 19, 2016

The point about cholesky was not very related to this PR. Instead it was related to the idea of blocked linear algebra in general. Cholesky is a very useful algorithm for common applications, probably much more useful than inv.

@mrocklin
Copy link
Member

mrocklin commented Jan 19, 2016

cc @jcrist and @ahmadia who might find this PR of interest

@ahmadia
Copy link

ahmadia commented Jan 19, 2016

This is a very good algorithm to test "out-of-core", since the memory complexity is going to stay at O(n^2), where n is the dimension of the matrix. This looks good to me, I wish I had dask when I was teaching parallel computing to HPC folks!

@sinhrks sinhrks force-pushed the solve_tri branch 2 times, most recently from cfd3792 to e91f4cf Compare January 21, 2016 18:39
@sinhrks
Copy link
Member Author

sinhrks commented Jan 21, 2016

Rebased on #885.

mrocklin
mrocklin reviewed Jan 21, 2016
View reviewed changes
dask/array/linalg.py Outdated
Copy link
Member

@mrocklin mrocklin Jan 21, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Need dtype here as well

Copy link
Member Author

@sinhrks sinhrks Jan 21, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've left because it is likely to be float, but not always. Isn't it the same as other funcs?

Copy link
Member

@mrocklin mrocklin Jan 21, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh, indeed. I should change that. In almost all cases we include the dtype. Some of the downstream libraries, like xray, find this very valuable. I'll fix the other linalg functions.

In test_array_core.py the eq function actually tests that dtype is set.

@mrocklin mrocklin merged commit 8bfc2a2 into dask:master Jan 22, 2016
@mrocklin
Copy link
Member

mrocklin commented Jan 22, 2016

This looks great. Merged.

@sinhrks sinhrks deleted the solve_tri branch January 22, 2016 19:19
@sinhrks sinhrks added this to the 0.8.0 milestone Feb 27, 2016
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Reviewers

No reviews

Assignees

No one assigned

Labels

None yet

Projects

None yet

Milestone

0.8.0

Development

Successfully merging this pull request may close these issues.

3 participants

@sinhrks @mrocklin @ahmadia