Is there a function/subroutine in MKL to specify that the product (matrix multiplication) of two different matrices would be symmetric so that matrix multiplication happens faster. I am aware of DSYR2K but that routine only helps when we are multiplying A*A' which would result in a symmetric matrix.
E.G. I want to multiply H*Q*H' where Q is positive definite and H is any arbitrary rectangular matrix. These types of matrix multiplications arise quite a lot in statistics and inverse problems. I know before hand that output of HQH' would be a symmetric matrix. Hence, when I multiply HQ= H*Q; B=HQ *H'; I want to tell MKL that B is a symmetric matrix to reduce the cost of matrix multiplication.
I understand that it is difficult for the software system to anticipate that the product of two matrices would be symmetric, but don't you think that it also applies to other routines. For ex in MKL_SCOOMM a user has to tell the routine that the left or /right hand side matrix would be symmetric (of course it has to be square but it need not be symmetric !) so that she/he can get some comupational benefit. Now in this case it is upto the user to make sure that the output would be a symmetric matrix.
In terms of above given terminology I can efficiently do H*Q from mkl but I cannot further reduce the cost of matrix multiplication by taking advantage of the fact that the output of HQ*H' would be a symmetric matrix. For small cases this is not even important, however it becomes an important issue when you are dealing with matrices that would be million by million in dimensions; of course such an operation would be done through Scalapack.
As far as the applications are concerned see paper (freely available online) http://www.geosci-model-dev-discuss.net/5/3325/2012/gmdd-5-3325-2012-print.pdf (the matrix multiplication algorithm is described for the Kronecker Product forms).For resons for dealing with these large dimensions see: http://onlinelibrary.wiley.com/doi/10.1029/2011WR011778/abstract ;
Lastly :The matrix multiplication of two sparse matrices has been raised earlier in the forum!: In some cases (~70 %) both (HQ) and H' are sparse matrices.
I understand that it is difficult for a software subsystem to know in advance that the output of the matrix multiplication of two matrices would result in a symmetric matrix. However, don’t you think so that a similar situation arises in case of using SSYMM as the routine does not know that the input left or right hand side matrix is a symmetric matrix (of course it would be square) and it is up to the user to make sure to input a symmetric matrix otherwise the output would be incorrect.
Now in the case mentioned by me in my earlier post it is also up to the user to make sure that the resulting matrix of the matrix multiplication would be a symmetric matrix.
A slight modification of DGEMM can fix this by only computing upper or lower triangular part of the matrix. Additionally, similar modification to sparse routines like MKL_DCSMM can fix this problem for sparse matrices. I have also raised this issue about BLAS in my recently submitted paper titled “Technical Note: Improving computational efficiency in large linear inverse problems: an example from carbon dioxide flux estimation” (Search on Google; freely available; see section 2.3).
I think, last time I did not answer your question in the context of real-time dynamic computing environment. I agree that there exists no method to find out in a real time dynamic computing environment if the output of the matrix multiplication of two arbitrary matrices would be a symmetric matrix. However, in case of three matrices where the central matrix is a symmetric matrix and left and right hand side matrices are transposes of each other this would always be true. For e.g. assume ‘a’ to be arbitrary and ‘q’ to be symmetric then a*q*a’ would be symmetric.
Let’s think of a simple application of computing regression coefficients where q is a full covariance matrix and a is a matrix of independent variables, then the most fundamental operation to get regression coefficient is to multiply a*q*a’ where q is mostly symmetric and these kinds of multiplications are quite common in statistics, econometrics, weather prediction, inverse modeling and maybe more‼
Many thanks for responding and initiating a healthy discussion
I completely agree with you that ~2.5% gain in performance is important especially when sometimes in my applications the matrices are huge reaching ~ million*million in dimensions in which case out of core operations becomes a necessity. I think the concept of
symmetry can be useful in improving performance in case of Strassen's algorithm and theoretically it is easy to show this performance gain. However, given the recursive nature of the algorithm programmatically designing it efficiently might be tricky.
For performance reasons in our applications, for toeplitz matrices we have now started using discrete fourier transform where we have to perform only o(n2logn) operations for obtaining the result of matrix multiplication.
Sometimes where absolute accuracy is not important we have also started using hierarchical matrices for matrix multiplication that again requires o(n2 logn) operations.
However both these methods can only be used in specific circumstances but still are extremely useful
(These algorithms I think are not available in MKL as they can only be used in very specific circumstances)