RFC: Matrix logarithm function (issue #5840) #12217

mfasi · 2015-07-19T17:41:21Z

I added the matrix logarithm function, with test cases. Although I tried to stick to the guidelines, it's my first contribution, please forgive any silly errors I might have committed.

joehuchette · 2015-07-19T17:54:41Z

base/linalg/dense.jl

@@ -278,9 +278,42 @@ function rcswap!{T<:Number}(i::Integer, j::Integer, X::StridedMatrix{T})
    end
 end

+function logm{T}(A::StridedMatrix{T})


The type parameter T isn't used anywhere?

It looks like it isn't, should I remove it from the signature of the function?

Unless you're using it 1) in the function body, or 2) to filter the types for which this method applies, then yeah, it's not necessary.

As mentioned by @joehuchette, the T parameter should be removed here.

jiahao · 2015-07-19T17:56:51Z

Ref #5840 - GitHub does not link issue numbers from titles

andreasnoack · 2015-07-19T18:21:32Z

Thanks for contributing this. It appears that (at least some of) the code is a translation of http://eprints.ma.man.ac.uk/1687/02/logm.zip. According to http://www.mathworks.com/matlabcentral/fileexchange/33393 it should be BSD licensed so there is no problem in using it, but I think you should link to the source in the initial comment.

andreasnoack · 2015-07-19T18:22:49Z

base/linalg/dense.jl

+    # Use Schur decomposition
+    n = chksquare(A)
+    S,Q,d = schur(complex(A))
+    R = full(logm(UpperTriangular(S)))


I don't think full is needed here.

tkelman · 2015-07-19T18:28:49Z

Note that there's the awful "Mathworks products only" clause in the File Exchange terms of use, see 2.a.iii. at http://www.mathworks.com/matlabcentral/termsofuse.html#content. So despite the BSD license on the code, unless it's available through a different medium or directly from the author (edit: with a clear license or permission to redistribute), having looked at any code hosted on File Exchange means we can't distribute a translation of it in Julia.

andreasnoack · 2015-07-19T18:36:11Z

The first link above is to one of the authors' homepage so hopefully Mathworks can't claim this, but I cannot find any license statement in that source.

andreasnoack · 2015-07-19T18:49:32Z

base/linalg/triangular.jl

+        A = sqrtm(A)
+    end
+
+    AmI = A - UpperTriangular(eye(n))


You can do this more efficiently with A - I. I is a special constant of type UniformScaling.

It looks like it doesn't work if A is UpperTriangular, but

UpperTriangular(full(A) - I)

is even slower.

Okay. Now I understand. This is supposed to work so I'll fix this now.

Thank you. It looks like scale(A,x) has the same issue when A is UpperTriangular.

ViralBShah · 2015-07-19T18:50:02Z

Perhaps we can write to the author and get clear documentation for the license.

andreasnoack · 2015-07-19T19:23:30Z

base/linalg/triangular.jl

+    end
+
+    AmI = A - UpperTriangular(eye(n))
+    d2 = norm(AmI^2, 1)^(1/2)


sqrt and cbrt are much faster

mfasi · 2015-07-19T19:41:59Z

I received the matlab code form the author, I think I can ask him for the permission to redistribute. Anyway, thank you for all your useful comments.

mfasi · 2015-07-20T13:00:59Z

The code is released under the terms of the BSD license. The author has just uploaded a new version of the .zip file containing the licence statement (http://eprints.ma.man.ac.uk/1851/02/logm.zip).

andreasnoack · 2015-07-20T15:10:11Z

base/linalg/triangular.jl

+        A = sqrtm(A)
+    end
+
+    AmI = UpperTriangular(full(A) - I)


This should be fixed now so you should be able to do AmI = A - I

andreasnoack · 2015-07-20T15:15:39Z

Great that the license question was resolved so fast. I've added a few comments, but I think it should be ready soon.

…ecking to ensure the precision to be double

andreasnoack · 2015-07-20T19:10:40Z

@mfasi Do you think this pr is ready now? The AppVeyor failure is unrelated.

mfasi · 2015-07-20T20:25:50Z

@andreasnoack All that I wanted to change is there and I'd be happy if you merged this pull request. Thank you for your feedback and support.

RFC: Matrix logarithm function (issue #5840)

andreasnoack · 2015-07-20T20:51:37Z

@mfasi Thank you for the contribution

tkelman · 2015-07-20T21:00:16Z

We do actually need to reproduce the original license statement in full here, for derived code.

Keno · 2015-07-20T21:03:08Z

We can put it in License.md with a reference to which code falls under the license and just a note of This code is derived from code (c) The Author 2015, see License.md for full details.

tkelman · 2015-07-20T21:08:53Z

There's a list of exceptions under LICENSE.md, and in contrib/add_license_to_files.jl, and I had modified a few of the license headers wherever files contained code derived from some other work. Or we could ask the author if they're willing to allow @mfasi to relicense his contribution under Julia's MIT license, then we could just put a comment indicating as much on this function.

tkelman · 2015-07-20T21:32:06Z

This is also resulting in warnings which should be suppressed

     * linalg3              WARNING: Matrix with nonpositive real eigenvalues, a nonprimary matrix logarithm will be returned. 
WARNING: Matrix with nonpositive real eigenvalues, a nonprimary matrix logarithm will be returned. 
WARNING: Matrix with nonpositive real eigenvalues, a nonprimary matrix logarithm will be returned. 
WARNING: Matrix with nonpositive real eigenvalues, a nonprimary matrix logarithm will be returned.

mfasi · 2015-07-20T22:07:47Z

@tkelman Sure, I can ask the author to relicense under Julia's MIT license. Is it enough if he just tells me that I can or should he add something to the original code?

I don't know how to suppress the warnings, nor I could find it anywhere. What should I do? The warnings are an expected behaviour in that case.

tkelman · 2015-07-20T22:11:53Z

If you get written permission, we can add a comment indicating "Julia version relicensed with permission from original author" and I think it should be fine (though IANAL).

Do the warnings occur every time, or only some of the time depending on randomness in the input matrices? There are other cases in the tests that do redirect_stderr, or we could set a seed so the input values are the same every time. It is good to test warning code paths for the sake of code coverage, but we should do so in a way that suppresses expected warnings wherever possible.

mfasi · 2015-07-20T22:22:37Z

I'll contact the author and let you know. There is no randomness involved, the warnings are produced every time you run the tests. Thus, I think that redirect_stderr would be fine in this case.

ViralBShah · 2015-07-26T11:26:09Z

@mfasi Thanks for patiently following the process and getting logm in. @higham Thank you for providing logm under the BSD/MIT licenses.

Cc: @alanedelman

tkelman · 2015-07-26T11:49:28Z

Just a hypothetical: how many tenured math, cs, and stats faculty around the world do you think we (+scipy, R, and octave communities) could get to co-sign a letter to the MathWorks complaining that the "may only be used with MathWorks products" clause is toxic (and hidden from users who otherwise think they're releasing code with a normal open source license) and asking them to drop it?

ScottPJones · 2015-07-26T11:57:19Z

Good idea. Do you have a link to the toxic clause in the MathWorks license?

mfasi · 2015-07-26T12:21:57Z

@ViralBShah You're very welcome. I'll add a link to this pull request in #5840.

joehuchette reviewed Jul 19, 2015
View reviewed changes

tkelman added the linear algebra Linear algebra label Jul 19, 2015

andreasnoack reviewed Jul 19, 2015
View reviewed changes

base/linalg/triangular.jl

end

AmI = A - UpperTriangular(eye(n))

d2 = norm(AmI^2, 1)^(1/2)

Copy link

Member

andreasnoack Jul 19, 2015

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

sqrt and cbrt are much faster

mfasi force-pushed the matrix_functions branch 2 times, most recently from 84090ed to d6add4a Compare July 20, 2015 10:21

Matrix logarithm function logm added

4d9bffa

andreasnoack reviewed Jul 20, 2015
View reviewed changes

mfasi force-pushed the matrix_functions branch from d6add4a to f056fda Compare July 20, 2015 15:14

mfasi force-pushed the matrix_functions branch from f056fda to 9763476 Compare July 20, 2015 15:29

A few minor improvements, high order functions removed, added type ch…

059d870

…ecking to ensure the precision to be double

mfasi force-pushed the matrix_functions branch from 9763476 to 059d870 Compare July 20, 2015 16:00

andreasnoack added a commit that referenced this pull request Jul 20, 2015

Merge pull request #12217 from mfasi/matrix_functions

7e7dfd3

RFC: Matrix logarithm function (issue #5840)

andreasnoack merged commit 7e7dfd3 into JuliaLang:master Jul 20, 2015

mfasi deleted the matrix_functions branch July 20, 2015 22:35

mfasi mentioned this pull request Jul 20, 2015

Fix warnings due to logm in tests/linalg3.jl #12236

Merged

tkelman mentioned this pull request Jul 24, 2015

ArgumentError: invalid argument #3 to LAPACK call in gebal! #12294

Closed

mfasi mentioned this pull request Jul 26, 2015

Implement more matrix functions #5840

Closed

17 tasks

argriffing mentioned this pull request Jul 31, 2015

fill all matrix function checkboxes in table 1 of http://eprints.ma.man.ac.uk/2102/ scipy/scipy#3365

Open

10 tasks

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

RFC: Matrix logarithm function (issue #5840) #12217

RFC: Matrix logarithm function (issue #5840) #12217

mfasi commented Jul 19, 2015

joehuchette Jul 19, 2015

mfasi Jul 19, 2015

joehuchette Jul 19, 2015

andreasnoack Jul 20, 2015

jiahao commented Jul 19, 2015

andreasnoack commented Jul 19, 2015

andreasnoack Jul 19, 2015

tkelman commented Jul 19, 2015

andreasnoack commented Jul 19, 2015

andreasnoack Jul 19, 2015

mfasi Jul 19, 2015

andreasnoack Jul 19, 2015

mfasi Jul 19, 2015

ViralBShah commented Jul 19, 2015

andreasnoack Jul 19, 2015

mfasi commented Jul 19, 2015

mfasi commented Jul 20, 2015

andreasnoack Jul 20, 2015

andreasnoack commented Jul 20, 2015

andreasnoack commented Jul 20, 2015

mfasi commented Jul 20, 2015

andreasnoack commented Jul 20, 2015

tkelman commented Jul 20, 2015

Keno commented Jul 20, 2015

tkelman commented Jul 20, 2015

tkelman commented Jul 20, 2015

mfasi commented Jul 20, 2015

tkelman commented Jul 20, 2015

mfasi commented Jul 20, 2015

ViralBShah commented Jul 26, 2015

tkelman commented Jul 26, 2015

ScottPJones commented Jul 26, 2015

mfasi commented Jul 26, 2015

RFC: Matrix logarithm function (issue #5840) #12217

RFC: Matrix logarithm function (issue #5840) #12217

Conversation

mfasi commented Jul 19, 2015

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jiahao commented Jul 19, 2015

andreasnoack commented Jul 19, 2015

Choose a reason for hiding this comment

tkelman commented Jul 19, 2015

andreasnoack commented Jul 19, 2015

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ViralBShah commented Jul 19, 2015

Choose a reason for hiding this comment

mfasi commented Jul 19, 2015

mfasi commented Jul 20, 2015

Choose a reason for hiding this comment

andreasnoack commented Jul 20, 2015

andreasnoack commented Jul 20, 2015

mfasi commented Jul 20, 2015

andreasnoack commented Jul 20, 2015

tkelman commented Jul 20, 2015

Keno commented Jul 20, 2015

tkelman commented Jul 20, 2015

tkelman commented Jul 20, 2015

mfasi commented Jul 20, 2015

tkelman commented Jul 20, 2015

mfasi commented Jul 20, 2015

ViralBShah commented Jul 26, 2015

tkelman commented Jul 26, 2015

ScottPJones commented Jul 26, 2015

mfasi commented Jul 26, 2015