Matrix norms implementation

[MpoMp]

Progress on #54! Read the issue comments for more info on the progress so far. The numbers in the tests were compared with the results of the same runs in MATLAB.

Looks like the singular values returned by MATLAB's svd and LAPACK's gesvd differ after some decimal places. Any feedback on that?

[MpoMp]

I converted the lambdas to protected functions. Now I am trying to manipulate yale and list matrices when calculating the one norm. Is there a specific method which returns the A array of a yale matrix? Seems like I haven't gotten the grasp of what is returned when calling all regular matrix functions of dense matrices.

[translunar]

Good. One more change, if you would. Best practice is to use full words (:infinity, :frobenius) for options instead of partials. I think it's okay to accept partials, but full words need to also work. After that I can merge. =)

Thanks!

[agarie]

Can you break this line into two? Keeping a limit of 80 characters per line is good practice!

[MpoMp]

Sure, I was wondering what's the exact limit anyway. Will commit this along with the planned changes I am working with now. :)

[MpoMp]

Done! I also added some TODO for @MohawkJohn's suggestions and some other points of the code which I am doubtful about.

[MpoMp]

Can I have some feedback please? :)
I think @MohawkJohn should be pleased with this. The only issues I can see now is making it work for sparse matrices and check if gesvd gives the desired values.

[agarie]

Pretty sure you can rearrange this giant if-else into a case expression:

case type
when nil then self.two_norm
when 1 then self.one_norm
when 2 then self.two_norm
when :frobenius, :fro then self.fro_norm
end

Play with it for some time. Ruby's case expressions are very cool. 🎉

[MpoMp]

I was thinking about this the other day. You're right, I'll push a commit in a few minutes. :)

[MpoMp]

Well, it looks sooo much better now. Thanks for the tip @agarie! :-)

[translunar]

Are you creating an NMatrix of arrays here? That's what it looks like, and it seems unsafe.

[MpoMp]

No, it's the NMatrix on which the norm calculation occurs, converted into a column vector. A column is appended per iteration. Should I change this to an 1xn NMatrix?

[MpoMp]

Fixed bad calculations and added some more. It was done based on SciPy's documentation for the relevant method: https://github.com/numpy/numpy/blob/v1.8.0/numpy/linalg/linalg.py#L1924

Should I reference their code somewhere? I didn't copy anything, I just tested my results against theirs. MATLAB totally avoided. :)

gesvd still returns a slightly different result (differs after some decimal places from the SciPy test results) from the desired one.

[translunar]

This is inefficient. Create the NMatrix but without default values, since you know its final shape, and insert the values directly. You're doing array construction, array resizing, and then NMatrix initialization from that array, which is an order of magnitude more expensive than necessary.

[translunar]

It's extremely bizarre that the results would differ from those in SciPy. What method are they using to perform this calculation?

And yes, I think if you do a #norm call on an NMatrix, it should return an NMatrix. That makes more sense than returning an Array.

[MpoMp]

Arrays exterminated! =)

Besides that, after doing some tests today, I discovered the svd results' bias might be a precision issue. Defaulting the self_cast matrix to :float64 gives a better result, though still not the same as SciPy's.
gesvd gives 7.3484692283495 35 while their test gives 7.3484692283495 45.

Furthermore, I guess there is an overflow issue when it comes to the lowest singular value because for :float64 it gives 0 instead of the really small number it should return (1.8628605857884395e-07 on my previous runs with :float32, 1.8570331885190563e-016 on the corresponding SciPy test). Any lights on this?

[abhinavagarwalla]

Is the feature fully implemented ? If it is not, I would really like to contribute.

[MpoMp]

I never got to finish it since I did not get any answer to my last questions and was unable to "dive deeper" on my own. Now I have other obligations and I guess I can't keep working on this. From my side: you're good to go. I guess @MohawkJohn will not have any objections either. :) Good luck!

[abhinavagarwalla]

@MohawkJohn Could you please guide me a little, because I am completely new to open source and is completely clueless at the moment on how to proceed. What all modifications can I do ?

[translunar]

@abhinavagarwalla Hey, sorry, somehow I missed your comment. Are you still working on this?

[yoongkang]

Is it really necessary to compare float64 results to 16th decimal place?

I think we should change the test to compare the difference between the expected result and and the computed result and accept it if it falls below a sensible threshold (say, 1.0e-10) and then call it a day.

Thoughts? I'd be happy to put that in if that's the only thing holding this up.

[yoongkang]

New PR at #400.

[v0dro]

Should this be closed? Outdated.