Skip to content
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

Sparse optimizations for GaBW sampling #331

Open
wants to merge 3 commits into
base: develop
Choose a base branch
from

Conversation

lucaperju
Copy link
Contributor

Samples faster when the A matrix of the polytope is sparse.
Very fast when both the A matrix of the polytope and the covariance matrix are sparse.

Copy link
Contributor

@vfisikop vfisikop left a comment

Choose a reason for hiding this comment

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

Great work thanks!

template <typename update_parameters>
auto compute_reflection(Point &v, Point &p, update_parameters const& params) const
-> std::enable_if_t<std::is_same_v<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>> && !std::is_same_v<update_parameters, int>, void> { // MT must be in RowMajor format
-> std::enable_if_t<std::is_same_v<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>> && !std::is_same_v<update_parameters, int>, void> {
Copy link
Contributor

Choose a reason for hiding this comment

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

What is update_parameters here? I do not understand !std::is_same_v<update_parameters could you please explain?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

basically there's another compute_reflection function above which takes an integer (just the facet) as the 3rd argument, and for some reason, if I don't have that condition the compiler decides to call this function assuming that the typename of update_parameters is int. There might maybe be better ways of dealing with these issues, but I remember I tried to solve them for some time and this is the best I could do, I couldn't at all understand how the compiler chooses which function to use when there's multiple ones that match

Copy link
Contributor

Choose a reason for hiding this comment

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

It seems that we are too generic here and that complicates the interface a lot! compute_reflection is called by several walks. Instead of creating all those complicated overloads why not simply create a new name. Especially if this is the case that this function is only called by a single walk (i.e. accelerated billiard walk). Moreover, update_parameters should not be a template but the struct defined in billiard walk, in all other cases this code will not compile since all other "update_parameters" does not have a moved_dist field.

include/convex_bodies/hpolytope.h Outdated Show resolved Hide resolved
include/convex_bodies/hpolytope.h Outdated Show resolved Hide resolved
typedef typename Polytope::VT VT;
typedef typename Point::FT NT;
using AA_type = std::conditional_t< std::is_same_v<MT, typename Eigen::SparseMatrix<NT, Eigen::RowMajor>>, typename Eigen::SparseMatrix<NT>, DenseMT >;
Copy link
Contributor

Choose a reason for hiding this comment

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

similar comment to AE_type, is there a better naming?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

hmm, I'll think of one, I'm not really sure what name I could give it but I'll see if I can come up with a better name.

if constexpr (std::is_same<AA_type, Eigen::SparseMatrix<NT>>::value) {
_AA = (P.get_mat() * P.get_mat().transpose());
} else {
_AA.noalias() = (DenseMT)(P.get_mat() * P.get_mat().transpose());
Copy link
Contributor

Choose a reason for hiding this comment

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

should we explicitly cast it to DenseMT?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I'm not sure what you mean, for the optimizations I need it to be in colmajor SparseMatrix format

Copy link
Contributor

Choose a reason for hiding this comment

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

the question is: "Where is the problem with _AA.noalias() = P.get_mat() * P.get_mat().transpose(); ?

E_type _E;
VT _AEA;
unsigned int _rho;
update_parameters _update_parameters;
typename Point::Coeff _lambdas;
typename Point::Coeff _Av;
bool was_reset;
BoundaryOracleHeap<NT> _distances_set;
Copy link
Contributor

Choose a reason for hiding this comment

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

This is defined in uniform ABW, it should be better if it is defined in a separate file and both walks include it.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

yeah, I was thinking about that too, any suggestions for that file name?

Copy link
Contributor

Choose a reason for hiding this comment

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

boundary_oracle_heap?

_distances_set.vec[i].first = ( *(b_data + i) - (*(Ar_data + i)) ) / (*(Av_data + i));
}
// rebuild the heap with the new values of (b - Ar) / Av
_distances_set.rebuild(_update_parameters.moved_dist);
Copy link
Contributor

Choose a reason for hiding this comment

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

Why not inserting the new values in the heap (in O(logn)) instead of rebuilding (in O(n))?

Copy link
Contributor Author

@lucaperju lucaperju Sep 13, 2024

Choose a reason for hiding this comment

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

here, this happens after we set a new direction, in which case now all the values have changed, so it's quicker to rebuild (O(n)) rather than insert each one (O(nlogn)). I'm not entirely sure if it makes a difference, but I think it does, since afterwards I never do O(nlogn) things, just O(non_zeroes * logn), so this O(nlogn) could be a bottleneck

@vfisikop vfisikop changed the title Spare optimizations for GaBW sampling Sparse optimizations for GaBW sampling Oct 9, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants