add relative_error and project_array functions#410
Conversation
| gramian = basis.gramian(product) | ||
| rhs = basis.inner(U, product) | ||
| coeffs = np.linalg.solve(gramian, rhs).T | ||
| return basis.lincomb(coeffs) |
There was a problem hiding this comment.
Is there a reason why not use basis = gram_schmidt(basis, product=product, atol=0, rtol=0) here? I would guess gramian could be ill-conditioned.
There was a problem hiding this comment.
Good point! A problem with gram_schmidt might be that one would have to make a copy of basis, which potentially requires a lot of memory. Not sure what to prefer ...
There was a problem hiding this comment.
The user could always call gram_schmidt outside and then set orthonormal to True. How about keeping the current code and additionally computing the condition of gramian? Then one could issue a warning, if the condition is too large. - On the other hand, the condition probably gets large very quickly. Any opinions, @pmli?
There was a problem hiding this comment.
@sdrave How about adding copy=False to the gram_schmidt call? (That would actually make more sense than my first suggestion.)
(BTW, since recently, scipy.linalg.solve raises a warning if the system is ill-conditioned.)
There was a problem hiding this comment.
copy=False is not really an option. The caller really would not expect basis to change, just by computing a projection.
There was a problem hiding this comment.
Aha, ok. Then replacing numpy.linalg.solve by scipy.linalg.solve would be good.
SciPy's version will warn about ill-conditioned gramians and should be slightly faster.
| else: | ||
| gramian = basis.gramian(product) | ||
| rhs = basis.inner(U, product) | ||
| coeffs = scipy.linalg.solve(gramian, rhs, sym_pos=True, overwrite_a=True, overwrite_b=True).T |
There was a problem hiding this comment.
sym_pos=True should probably be replaced by assume_a='pos' (see here). But actually, is gramian guaranteed to be numerically symmetric and positive definite? Adding gramian = (gramian + gramian.T) / 2 might be a good idea, and solve will probably detect if gramian is close to being singular.
There was a problem hiding this comment.
I used sym_pos instead of assume_a since the latter is only present in the latest SciPy release and I didn't want to make the code more complicated. - I guess in presence of a product, there is no guarantee that gramian is exactly spd. However, why do you think, that gramian = (gramian + gramian.T) / 2 would improve the result of the solve call in that case?
There was a problem hiding this comment.
I guess it's mostly for my own sanity. But, it's nice that (A+A^T)/2 is the closest symmetric matrix to A. Anyway, Cholesky factorization can be computed in a stable way, so there shouldn't be a huge difference. Probably the bigger issue is that the condition number of gramian is the square of the condition number of basis (at least in the standard inner product).
Anyhow, calling gram_schmidt before avoids all of this, so feel free to ignore my previous comments.
There was a problem hiding this comment.
Ok, I would also doubt that A+A^T/2 has an important effect. Merging as is ..
2 participants
No description provided.