Skip to content

WIP: Generalized matrix product and a fix in endomorphic MatrixProductOperator

Neel Shah requested to merge g-neelshah/nifty:generalized_matrix_product into NIFTy_8

I added a linear matrix multiplication operator that generalizes the endomorphic MatrixProductOperator. The output matches that of the endomorphic operator in case of an endomorphic application.

The target domain/domain tuple of the valid shape should be specified by the user, but I kept a standard RGSpace of the valid shape as the default.
Below is an example of an application.

dtuple = DomainTuple.make(((RGSpace(10),RGSpace(20),RGSpace(30),RGSpace(40)))
spaces = (1,3) (apply only to 2nd and 4th spaces of domain)
matrix = np.ones(shape=(5,15,25,20,40))
matop = GeneralMatrixProduct(dtuple,matrix,spaces,target=None)

According to the operator's convention,
target_shape = (10,5,30,15,25)
(since dtuple.shape = (10,20,30,40), the 20,40 axes are summed over, and to comply with the endomorphic MatrixProductOperator the 10,30 axes retain their positions and remaining axes of the matrix are filled in order) = DomainTuple.make(RGSpace(shape = target_shape, distances = default, harmonic = False))
(unless a domain or DomainTuple with shape = target_shape is passed as the target)
returns Field.from_raw(, np,tensordot of matrix and field summed over appropriate axes)

I tested the operator's working on the example above, and that it coincides with the endomorphic MatrixProductOperator when the latter is applicable. Would this be a useful operator, and what else in it can be tested or changed?

The second change is in the endomorphic MatrixProductOperator, I fixed its functionality for application to a >1 dimensional domain (say an RGSpace(shape=(5,15,25))). The previous version uses instead of np.tensordot, so due to default axis convention in it doesn't work when more than one axes are to be summed over. I don't think this was intentional because the operator doesn't throw an error until it is applied to a field over such a domain.

There is one last change in the docstring of the matrix product operator. It has a line saying If DomainTuple it is assumed to have only one entry, which I think is misleading because I can apply it to a domain tuple with many subspaces like the example above. I changed the line to If Domain it is assumed to have only one subspace, which is what I observed.

@gedenhof , can you have a look at this?

Merge request reports