For two matrices A and B the trace(A^T*B) is just Aij*Bij with implied summation. Is (A.transpose()*B).trace() the right way to code this? A.transpose().lazyProduct(B).trace()?

Or is there some better expression based on the result just being a dot product of the coefficients when viewed as vectors? E.g. for two Matrix3d A & B, short of using a Map or {A.B).data(), is there a crisp way to write that dot product? A crisp way that accounts for possible differences in storage order?

Thanks,
Rhys

Possibly helpful, this operation is also called the Frobenius inner product (http://en.wikipedia.org/wiki/Matrix_mul ... er_product).

What I wanted was A.cwiseProduct(B).sum() or (A.array()*B.array()).sum(). Any better suggestions?

I don't have better suggestion, I can only add that among the 4 versions you wrote, only one is not recommended at the moment:

(A.transpose()*B).trace()

because A.transpose()*B will be fully evaluated into a temporary.

As a follow up I have a related question, are these the recommended construction for the corresponding linear algebra

A) (x'b-b0)^2
int n=1000;p=20;
MatrixXf x = MatrixXf::Random(n,p);
VectorXf b = VectorXf::Random(p);
float basecoef=0.263;
VectorXf proj = ((x*b).array()-basecoef).square();

B) (X-x_2_m)'S^{-1}_2(X-x_2_m)

X_2 and X are matrices,
x_2_m is the colwise mean of a matrix x_2
and S_2^{-1} is the inverse of the covariance matrix of x_2

int n=1000;p=20;
MatrixXf x = MatrixXf::Random(n,p);
MatrixXf x_2 = MatrixXf::Random(n,p);
RowVectorXf x_2_mean=x_2.colwise().mean();
x.rowwise() -= x_mean;
MatrixXf b = MatrixXf::Identity(p,p);
((x_2.transpose()*x_2)/(n-1)).ldlt().solveInPlace(b);
return ((x*b).array()*x.array()).rowwise().sum();

thanks in advance,

I don't understand what _2 means in B), I guess it's equivalent to:

trace( A' inv(S) A)

??

so I'd avoid to explicitly compute the inverse:

MatrixXf b = ((x_2.transpose()*x_2)/(n-1)).ldlt().solve(x);
return (x.array() * b.array()).rowwise().sum();

ggael wrote:I don't understand what _2 means in B), I guess it's equivalent to:

trace( A' inv(S) A)

??

so I'd avoid to explicitly compute the inverse:

MatrixXf b = ((x_2.transpose()*x_2)/(n-1)).ldlt().solve(x);
return (x.array() * b.array()).rowwise().sum();

for B I'd avoid to

Dear ggael,

thanks for your advice for B): i did not know that my construction was explicitly computing the inverse. (i have the impression part of your answer may have been lost in a copy paste).

I'm not sure this will work:

MatrixXf b = ((x_2.transpose()*x_2)/(n-1)).ldlt().solve(x);

((x_2.transpose()*x_2)/(n-1)) is a p by p matrix whereas x is n by p.
At any rate, the compiler throws "Assertion `m_matrix.rows()==b.rows() && "LDLT::solve(): invalid number of rows of the right hand side matrix b"

For the background i should have given before: X_2 is another matrix with the same .rows() as X. In the snipped of code above, i'm computing the Mahalanobis distances of the rows of X wrt to the mean/covariance of X_2 (a comon operation in statistical applications). From a linear algebra point of view it is not trace( A' inv(S) A) [A' inv(S) A is not square]. It is (A'inv(S)A).rowwise().sum().

sorry, I read too fast, you actually want the diagonal of x inv(S) x', so

ArrayXf b = ((x_2.transpose()*x_2)/(n-1)).ldlt().solve(x.transpose());
return (x.array() * b.transpose()).rowwise().sum();

or:

MatrixXf b = ((x_2.transpose()*x_2)/(n-1)).ldlt().solve(x.transpose());
return (x.lazyProduct(b)).diagonal();

Dear ggael,

i first turned back to my original solution,
but now i would like to understand why the
solution you proposed doesn't work, i.e.:


int h=100,n=200,p=5;
MatrixXf x_1=MatrixXf::Random(n,p);
MatrixXf x_2=MatrixXf::Random(h,p);
RowVectorXf x_2_mean=x_2.colwise().mean();
x_1.rowwise()-=x_2_mean;
x_2.rowwise()-=x_2_mean;
ArrayXf b=((x_2.transpose()*x_2)).ldlt().solve(x_1.transpose())*(h-1);
return (x_1.array()*b.transpose()).rowwise().sum();


returns:


PlainObjectBase.h:300: void Eigen::PlainObjectBase<Derived>::resizeLike(const Eigen::EigenBase<OtherDerived>&) [with OtherDerived = Eigen::CwiseUnaryOp<Eigen::internal::scalar_multiple_op<float>, const Eigen::ReturnByValue<Eigen::internal::solve_retval_base<Eigen::LDLT<Eigen::Matrix<float, -0x00000000000000001, -0x00000000000000001>, 1>, Eigen::Transpose<Eigen::Matrix<float, -0x00000000000000001, -0x00000000000000001> > > > >, Derived = Eigen::Array<float, -0x00000000000000001, 1>]: Assertion `other.rows() == 1 || other.cols() == 1' failed.
Aborted


In advance, thanks for your comments,

its:

ArrayXXf, not ArrayXf

thanks! (also for your patience)