Norm method format

[minivan]

Hey everyone,

I'd like to implement a norm method for NMatrix and would like to ask for your assistance.
How do you see the function? At the moment, here's what I think:

mat = NMatrix.new(4, 4)
mat.norm       # norm with p = 2
mat.norm 3     # norm with p = 3
mat.norm :inf  # infinity norm
mat.norm :fro  # Frobenius norm

What do you think?

[agarie]

@minivan I know little about matrix norms, but I think that the suggested API is interesting. A few questions:

• How can it be extended, i.e. how to add new norms?
• Can you give some use cases for matrix norms? This way, preparing the tests will be easier. Also, we'll know what NMatrix can do after it's implemented.

Go ahead. When you have some code or doubt about NMatrix's API, please send a message to the mailing list.

Thank you!

[translunar]

I like this proposal. I'd say :infinity and :frobenius for the symbols -- better to use the full names. Maybe also allow it to accept :inf and :fro.

[agarie]

Hey @minivan, are you still interested in this? :)

[minivan]

@agarie absolutely, but I got a bit stuck with porting a LAPACK function and getting the data types all correct for the dlange routine. I assume I don't have the proper libraries set up correctly, but I'll check everything now and see what help I can ask for! Thanks guys for your interest, I'm jumping right back onto it 🤘

[lazywei]

@minivan
Hi, I would like to get involved into this issue.
Could you please let me know which branch of your fork is about this issue?
Thanks!

[agarie]

@lazywei Looking at his fork, I don't see a branch that wasn't created by me or by John. So, I think you could start something on your own* and post your progress here (or as a pull request).

• Of course, feel free to ask things in IRC (#sciruby on freenode) or on the mailing list. :)

[translunar]

No work has been done on this in six months. Is anyone interested in finishing this? Otherwise I'm closing.

[MpoMp]

I am interested! I agree with the format proposed by the initial post:

mat = NMatrix.new([2, 2], [1, 2, 4, 8])
mat.norm       # norm with p = 2
mat.norm 3     # norm with p = 3
mat.norm :inf  # infinity norm
mat.norm :fro  # Frobenius norm

Symbols for norm names would be also accepted by their full names and case insensitive. I intend to implement the p-norm, frobenius (which would run the same code for p-norms) and max norm. Is that sufficient for a first stage?
Should I post in the mailing list or here when I have made progress? I think my code will not conform 100% to the guidelines and I am really in the dark when it comes to code testing in ruby (reading about this right now though).

[translunar]

As an FYI, the 2-norm has been implemented for vectors:

nmatrix/lib/nmatrix/math.rb

Line 535 in 18743e0

def nrm2 incx=1, n=nil

It's partially tested here:

nmatrix/spec/blas_spec.rb

Line 140 in 55417c3

it "exposes nrm2" do

but note that this is not the way we want future tests to be written. Those should look more like this one:

https://github.com/SciRuby/nmatrix/blob/master/spec/00_nmatrix_spec.rb#L476

(and note also that as I linked that I discovered the one after that is not properly named).

I think your proposal sounds good. Go for it. =)

[MpoMp]

My first attempt has been pushed! I added my implementation at the math.rb file, thought it was better than writing to the core NMatrix file. Check my nmatrix fork for the whole thing. Here are the gist pastes:
https://gist.github.com/MpoMp/9433134

I compared my output with the output at MATLAB. I might have mistaken something in the way I implemented. In all cases of the infinity norm my output was the same with MATLAB's. But when it comes to p-norms the results differ after the second decimal place. So I had to "hard-code" the values at the relevant tests.

So, is MATLAB a good way to test the output? I couldn't find a way to test p-norms for values over 2. Seems like I miss something from the relevant theory because they are not commonly implemented.

[translunar]

Couple of thoughts:

• Tie it in with the existing :norm2 method.
• After you get it working consistently in Ruby, try implementing it in C. It's going to be too slow in pure Ruby.
• It's very odd that you'd get a different answer from Matlab. Try GSL, perhaps, or R, or something else?
• Instead of sum**(1/p.to_f), try sum**(1.quo(p)).
• type.is_a?(Fixnum) is better than type.class == Fixnum because the former covers types derived from Fixnums as well, and is also implemented in C. Same with the other type.class comparisons.
• Converting symbols to strings is expensive, particularly when you then run a comparison. The advantage to sym == :blah is that it takes about 1 operation, whereas str == 'blah' may require four or more operations in the worst case (one for each character) — not to mention the cost of the conversion itself.
• You can use hashes to call functions or lambdas based on symbols more efficiently, I believe.
• Not a bad use of lambda, but these should probably be protected functions instead.
• In the future, I suggest using the pull-request mechanism instead of a gist. =)

In general, a very good first attempt! Hopefully these tips are useful.

[MpoMp]

1. Worked on using nrm2 within norm today. I doubt if it can be done with my current implementation because when a matrix norm is computed, the root calculation comes after both sums (inner and outer) have been calculated. If I use methods which return vector norms, the root calculation will be already applied leading to undesired numbers.
2. That would be implemented within ruby_nmatrix.c? And, if I get it right, the ruby part will just be calling the C function? I only have basic skills in C so unless you think that after reading through some guidelines I'll be able to give it a shot, I doubt I am the right person to attempt this!
3. Well, I had to read R's manual to discover my mistake. The 2-norm (when it comes to p-norm for p = 2) is not the same thing with the Frobenius norm. The p-norm is only applied to vectors. Neither MATLAB nor R implement it for matrices. I will change the norm method format to one similar to R's:

mat = NMatrix.new([2, 2], [1, 2, 4, 8])
mat.norm       # the “spectral” or 2-norm, which is the largest singular value of mat
mat.norm 2     # same as above
mat.norm 1     # one norm (absolute column sum)
mat.norm :inf  # infinity norm (absolute row sum)
mat.norm :fro  # Frobenius norm (the Euclidean norm of mat treated as if it were a vector)

The implementation of vector norms should be independent from this.
4. Done! It doesn't affect the results so I suppose it's a performance boost?
5. Done! Along with strings converted to symbols via hashes.
6. Protected functions outside the norm method?

Towards which branch should I request the pull? I'll request it on the master branch for now. I used gists because I didn't feel my work was worthy enough to appear in the main repository before you reviewed it. But I guess you can just reject the request if that's the case?

Thanks for the feedback! :)

[translunar]

Yes. Protected functions outside of the norm method.

Another concern I have is that you're doing column traversal for the 1-norm. This will be a fairly efficient operation for dense, but inefficient for list and especially yale. It may be that you wish to simply copy-transpose yale and list matrices and calculate the norm on the basis of rows instead.

[MpoMp]

Cool, I'll work on these two things as soon as possible and drop a line at the pull request when I have news.