Numerical sensitivity of function kron(X,Y)
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Assume X and Y are both symmetric matrices with double floating-point format. Can kron(X,Y) potentially become asymmetric due to the floating point numerics similar to the situation like X'*Y*X ?
I personally think it would be quite unlikely to have the problem I stated above given the mathematical definition of the Kronecker product. Nevertheless, if anyone has suggestions, please share them.
Thanks
0 commentaires
Réponse acceptée
Matt J
le 19 Déc 2016
Modifié(e) : Matt J
le 19 Déc 2016
My instinct and my tests say no. There are only scalar multiplications in A=kron(X,Y), so those should be pretty safe. However, if you plan to start multiplying A with other data, the benefits of the symmetry might go out the window pretty fast. Note, for example, that C as computed with,
C=reshape( kron(X,X)*Y(:) ,N,N);
should be symmetric in exact arithmetic for symmetric X and Y. However, I find that this fails in most numerical tests, likely due to the very same issue you pointed out with X.'*Y*X.
0 commentaires
Plus de réponses (1)
John D'Errico
le 19 Déc 2016
Never trust the LSBs in a floating point computation. If it can go bad, it probably will.
A = randn(200);
B = randn(100);
A = A + A';
B = B + B';
C = kron(A,B);
all(all(C == C'))
all(all(C == C'))
ans =
logical
1
I did my best above, but symmetry remains. C is pretty large, enough so that the blas should be kicking in, just in case they might have gotten me in trouble. I even chose different sizes for A and B, in case that might get me in trouble.
So I did a quick look at the guts of kron. It looks like it is safe to retain symmetry to me. A computation like X'*Y*X is a problem, because first MATLAB forms X'*Y, which is not symmetric, even if X and Y are symmetric. Then, even though a final multiplication by X should restore symmetry, that simply won't happen in floating point arithmetic.
Though still a conjecture and not proved, I think kron is ok.
2 commentaires
Voir également
Catégories
En savoir plus sur Logical dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!