grzegorzmazur / yacas

Computer calculations made easy
http://www.yacas.org
GNU Lesser General Public License v2.1
126 stars 24 forks source link

Symbolic IsSymmetric #327

Closed jeksterslab closed 3 years ago

jeksterslab commented 3 years ago

In the example below, I am symbolically deriving several symmetric matrices, namely Sigma which is a covariance matrix, Diag which is the inverse of the diagonal matrix consisting of the diagonal elements of Sigma, and R which is a correlation matrix. Shouldn't IsSymmetric return True for all matrices? The arrangement of variables may be different for Sigma and R but the lower and upper triangular elements are symbolically equivalent.

Input

A := {{0,alpha,1},{0,0,0},{0,0,0}};
B := {{0,0,0},{0,beta,0},{0,0,gamma}};
C := {{1,0,0},{0,1,0}};
Iden := Identity(3);
Sigma := C * Inverse(Iden - A) * B * Transpose(Inverse(Iden - A)) * Transpose(C);
Diag := Inverse(DiagonalMatrix(Sqrt(Diagonal(Sigma))));
R := Diag * Sigma * Diag;
IsSymmetric(B);
IsSymmetric(Diag);
IsSymmetric(Sigma);
IsSymmetric(R);

Output

{{0,alpha,1},{0,0,0},{0,0,0}}
{{0,0,0},{0,beta,0},{0,0,gamma}}
{{1,0,0},{0,1,0}}
{{1,0,0},{0,1,0},{0,0,1}}
{{beta*alpha^2+gamma,alpha*beta},{beta*alpha,beta}}
{{1/Sqrt(beta*alpha^2+gamma),0},{0,1/Sqrt(beta)}}
{{1,(alpha*beta)/(Sqrt(beta)*Sqrt(beta*alpha^2+gamma))},{(beta*alpha)/(Sqrt(beta*alpha^2+gamma)*Sqrt(beta)),1}}
True
True
False
False
grzegorzmazur commented 3 years ago

Hi,

Yeah, you've stumbled upon a serious issue in the very core of yacas, which is not transforming the expressions to the well-defined canonical form. As a result, in

{{beta*alpha^2+gamma,alpha*beta},{beta*alpha,beta}}

the off-diagonal terms have different order of alpha and beta. Hence they are not considered equal.

As a quick workaround I'd suggest using

Simplify(Transpose(A)-A) = 0

to check A symmetry. It involves Simplify() which should deal with the ordering issues just fine, but which in general case can be very computationally expensive - that's why I'd rather not put it in the standard IsSymmetric() predicate.

Cheers, Grzesiek

jeksterslab commented 3 years ago

Thanks very much.