Hi all,

I'm using Eigen and MATLAB to compute the norm of a matrix, and they're giving me slightly different answers. I believe MATLAB computes the norm of a matrix P by taking the square root of the largest eigenvalue of P'P (P' = transpose of P). However, Eigen must be doing differently. The following code

Code: Select all

   Matrix3f P;
   P << 1, 2, 3, 4, 5, 6, 7, 8, 9;
   Matrix3f Ptrans = P.transpose();
   Matrix3f PtransP = Ptrans*P;

   std::cout << "P'P = " << std::endl << PtransP << "\n\n";
   std::cout << "eigvals of P'P are " << std::endl << PtransP.eigenvalues() << "\n\n";
   std::cout << "sqrt(283.859) = " << sqrt(283.859) << "\n\n";
   std::cout << "norm of P is " << std::endl << P.norm() << "\n\n";

gives output

Code: Select all

P'P =
 66  78  90
 78  93 108
 90 108 126

eigvals of P'P are
(283.859,0)
(1.14142,0)
(-3.83749e-006,0)

sqrt(283.859) = 16.8481

norm of P is
16.8819


Clearly 16.8481 != 16.8819. Can anyone tell me which is the more accurate answer, and why?

Eigen provides stableNorm and blueNorm. These are supposed to increase accuracy/precision at the expense of speed. You might see if these give you a closer answer than norm.

graphicsMan: Good thought, but unfortunately I got this error when I tried to invoke stableNorm() on a matrix:

c:\program files\eigen\eigen-eigen-599fc4f7c736\eigen\src\core\vectorblock.h(135): error C2338: YOU_TRIED_CALLING_A_VECTOR_METHOD_ON_A_MATRIX

I guess stableNorm() only applies to vectors.

I seem to keep answering my own questions. The correct method was operatorNorm().

yes, and regarding 16.8481 != 16.8819 this is probably because you used float in Eigen while MatLab uses double.

@ggael: Actually was using doubles in Eigen too. My best guess is norm() uses some faster, but slightly less accurate, algorithm for computing the norm of a matrix. operatorNorm() gets it right on, but it's slower.

norm() computes the Frobenius norm: the square root of the sum of all the matrix entries. For your matrix, P.norm() = sqrt(1 + 4 + 9 + 16 + ... + 81).

@jitseniesen: I find that odd, considering the documentation states that norm() is the l2-norm of the matrix, which I interpret to be the spectral norm of a matrix (not simply the Euclidean norm of the matrix as a vector in Rmn). If you're right, I'd say the documentation should be changed (unless I am totally off-base in thinking that).

you are right, the documentation has to be more explicit.