-
Notifications
You must be signed in to change notification settings - Fork 12
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.
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
Facilitate streaming (row-by-row update) in ESMDA #162
Conversation
tests/test_esmda_inversion.py
Outdated
ans = function(alpha=alpha, C_D=C_D, D=D, Y=Y, X=X) | ||
K = function(alpha=alpha, C_D=C_D, D=D, Y=Y, X=None, return_K=True) | ||
|
||
X - np.mean(X, axis=1, keepdims=True) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Remove this
X - np.mean(X, axis=1, keepdims=True) | |
X - np.mean(X, axis=1, keepdims=True) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
LGTM.
This is very nice! In particular, see comparison between low-level and high-level API in tests. Good job!
Note for later: the names inversion
are slightly missleading, they all involve inversion, but their purpose is either to return K
or multiply X
with K
. Perhaps we can come up with good naming. As discussed, they are not part of the public api, so we may think of it later.
To support row-by-row (parameter by parameter) computation, either literally by updating a single parameter or updating batches, we should avoid doing duplicate work.
Consider the equation:
If we update row by row and introduce a transition matrix$K := Y^T /(N-1) (C_{DD} + \alpha C_D)^{-1} (D - Y) $ ,$[x_1 | x_2 | \ldots]^T K = [x_1 K | x_2 K | \ldots ]^T$ . In other words, we can pre-compute $K$ , since it's independent on $X$ , and apply it row-by-row or group-by-group, instead of re-computing it for each group.
then
This PR proposes two methods for this:
compute_transition_matrix
, which computes the transition matrixperturb_observations
, called by bothcompute_transition_matrix
andassimilate
, to avoid duplicating codeNote: this PR assumes #161 will be merged.
I also realized that we don't have to center both$X$ and $Y$ ! So I changed the code away from that, saving both memory and time by only centering the smaller matrix $Y$ . See this comment. This is also documented in the code with a few lines.