uncomplicate / neanderthal

Fast Clojure Matrix Library
http://neanderthal.uncomplicate.org
Eclipse Public License 1.0
1.06k stars 56 forks source link

issue with ev! when eigenvalues have multiplicity #49

Closed paroda closed 5 years ago

paroda commented 5 years ago

I was trying the examples from the tutorial: https://dragan.rocks/articles/17/Clojure-Linear-Algebra-Refresher-Eigenvalues-and-Eigenvectors

My env: [uncomplicate/neanderthal "0.20.4"] OS: Debian 9 MKL: 2019 Initial release

Noticed a few differences:

  1. The signature of ev! is now different from what is shown in example. Now it takes 4 arguments while in example it is shown taking 3 args.
  2. The calculation seems to be erroneous. I tried the same but getting different results. Also, the last test -Ax+λb = 0 is failing as can be seen below: (however, it works fine if tested on an example which doesn't have multiplicity > 1)
(require '[uncomplicate.neanderthal
           [core :refer [col entry mv scal axpy]]
           [native :refer [dge]]
           [linalg :refer [ev!]]])

(let [a (dge 3 3 [5 4 2, 4 5 2, 2 2 2])
      vr (dge 3 3)
      w (dge 3 2)]
  (ev! a w nil vr)
  [(col w 0)
   vr
   (let [l (entry w 0 0)
         v (col vr 0)]
     (axpy -1 (mv a v) (scal l v)))])

;; output =>
;;
;; [#RealBlockVector[double, n:3, offset: 0, stride:1]
;; [   1.00   10.00    1.00 ]
;;  #RealGEMatrix[double, mxn:3x3, layout:column, offset:0]
;;    ▥       ↓       ↓       ↓       ┓    
;;    →      -0.75    0.67    0.00         
;;    →       0.60    0.67   -0.45         
;;    →       0.30    0.33    0.89         
;;    ┗                               ┛    
;;  #RealBlockVector[double, n:3, offset: 0, stride:1]
;; [   0.00   -5.37    0.00 ]
;; ]
blueberry commented 5 years ago

The API has changed slightly since the time I write the blog post. If you want the examples to work 100% the same as in the text, it is best to use the version of neanderthal that was current at the time of writing the blog post, since I can't update all blog posts when there is a breaking change.

However, such breaking changes are not that numerous, and there is an extensive test suite and (I hope) sufficient documentation to help readers sort this out.

Now, related to this particular example, this is the direct copy-paste from my repl:

(let [a (dge 3 3 [5 4 2, 4 5 2, 2 2 2])
      vr (dge 3 3)
      w (dge 3 2)]
  (ev! a w nil vr)
  [(col w 0)
   vr
   (let [l (entry w 0 0)
         v (col vr 0)]
     (axpy -1 (mv a v) (scal l v)))])
[#RealBlockVector[double, n:3, offset: 0, stride:1]
[   1.00   10.00    1.00 ]
 #RealGEMatrix[double, mxn:3x3, layout:column, offset:0]
   ▥       ↓       ↓       ↓       ┓    
   →      -0.75    0.67   -0.03         
   →       0.60    0.67   -0.42         
   →       0.30    0.33    0.91         
   ┗                               ┛    
 #RealBlockVector[double, n:3, offset: 0, stride:1]
[   0.00   -5.37   -0.00 ]
]

You'll notice that it the result is the same as in the blog post.

Now, for the -5.37 part, that is, checking whether this is really the right result, I don't have time for this now, since I have to go outside, and, since I am not a matematician working with eigenvectors every day, I can't give you an answer from the top of my head, and would have to think about it and perhaps look things up.

Anyway, the heavy lifting is done by Intel MKL, and I doubt there is a bug there that would manifest with such a simple example. The key is in interpreting the results having in mind math theory...

paroda commented 5 years ago

Thanks for replying. I agree about the blog updates. It was minor. I am currently trying to learn linear algebra for deep learning and your library looks the most promising of all that i have tried so far.

And I do confess, I do not know much of the math. just starting. Should read more on the the theory about the case where eigenvalues have multiplicity.

For now, my main concern is to ensure that i have the right setup and so that i can proceed further, hence was confused with the results. Now, I tried again on another machine (win 7), and this time I got the same results as you have shown. So probably i got the setup right this time, and looks like i should be good to go.

blueberry commented 5 years ago

In this example, the problem you're running into with checking the parameters is that ev! is destructive. It overwrites a with the eigenvectors, so your test does not do what you intend, although the computed eigenvectors and eigenvalues are correct.

You need to preserve the original value of the matrix a, either by calling ev! on the copy of a, or by using the ev function.