PolyMathOrg / PolyMath

Scientific Computing with Pharo
MIT License
170 stars 39 forks source link

Symmetric matrix is not symmetric #37

Open olekscode opened 6 years ago

olekscode commented 6 years ago

We have two ways of instantiating PMSymmetricMatrix:

  1. by using the rows: method, inherited from PMMatrix
m := PMSymmetricMatrix rows: #(
    (1 2)
    (2 3)).
  1. by converting a PMMatrix to PMSymmetricMatrix using asSymmetricMatrix.
m := PMMatrix rows: #(
    (1 2)
    (2 3)).

m := m asSymmetricMatrix.

Neither of these methods checks whether a given matrix is really symmetric. This results in a better performance, especially in case of large matrices, when doing such test every time we create (or convert) a matrix could be expensive.

However, PMMatrix implements a method isSymmetric

isSymmetric
    ^ self = self transpose 

which is overridden by PMSymmetricMatrix to always return true

isSymmetric
    "Answers true because the receiver is a symmetric matrix"
    ^true

Which gives some confusing results

m := PMSymmetricMatrix rows: #(
    (1 0)
    (2 3)).

m isSymmetric. "true"
m := PMMatrix rows: #(
    (1 0)
    (2 3)).

m := m asSymmetricMatrix.
m isSymmetric. "true"

I think that either we should not allow initializingPMSymmetricMatrix with non-symmetric matrix (which is confusing already), or if we do (for the sake of performance), we should at least leave the method isSymmetric for testing if a symmetric matrix is actually symmetric.

nicolas-cellier-aka-nice commented 6 years ago

Hi Oleg, I think that the test should be performed at creation only, if ever it needs to be performed. Maybe it will require some tolerance, because for example, there is no guaranty that A*A transpose be numerically symmetric depending on the order of operations, and even less so for A*B transpose + (B*A transpose), but that's another question...

If we guaranty the symmetry by construction, then in no way should we modify the isSymmetric test. That would be only a waste of time and efficiency.

Other consideration might be to store only half the contents, either packed or not, like in LAPACK. There would be no way to test for symmetry after construction. So for me, I would not change the test.

olekscode commented 6 years ago

Hello, @nicolas-cellier-aka-nice , So should I add a test on creation? Because right now it's not performed

nicolas-cellier-aka-nice commented 6 years ago

IMO, we should have two creation messages, one with the test, the other without. In the later case it's user responsibility to ensure that the property is fulfilled. Generally, the matrix will be symmetric (or it's a programming bug), so the cost of checking will be n*(n-1)/2. In most cases, I think that the matrix will be symmetric by construction and the test will be superfluous. Another option would be to have some "debug" option that would assert some precondition here and there when true...

olekscode commented 6 years ago

eigen() function in R has a symmetric parameter

if TRUE, the matrix is assumed to be symmetric (or Hermitian if complex) and only its lower triangle (diagonal included) is used. If symmetric is not specified, isSymmetric(x) is used.

So it's similar to what you suggest with two methods:

matrix asSymmetricMatrix. "test is performed"
matrix asSymmetricMatrixNoTest. "test is not performed"