New LANumerics package for numerical linear algebra in Swift

Just open-sourced LANumerics, which is a Swift package for numerical linear algebra built on top of Swift-Numerics: GitHub - phlegmaticprogrammer/LANumerics: LANumerics is a Swift package for doing numerical linear algebra. .

It was no trouble at all to use Swift-Numerics! It would be good though to have better printing for complex numbers in Swift-Numerics, although I can work around that when they appear inside of matrices.

Also, it might make sense to change the definition of magnitude for complex numbers from max(real.magnitude, imaginary.magnitude) to real.magnitude + imaginary.magnitude as that is what is supported by the BLAS function for finding the index of the largest element.

8 Likes

Looks nice.

Hereā€™s some feedback:

First, have you considered making Element unconstrained, and putting protocol-specific operations in constrained extensions of Matrix? I personally find that yields a cleaner implementation and interface.

Have you considered making adjoint and transpose effectively lazy? As in, make new types: AdjointMatrix and TransposeMatrix, so that the elements do not actually get modified nor rearranged in memory.

I notice that youā€™ve included BLAS, LAPACK, and VDSP functions with un-Swifty names and signatures. Have you considered wrapping those in simpler, more natural spellings?

Speaking of naming conventions:

Swift uses description, not toString().
Swift uses initializers for type conversions, eg. String(precision:)
Swift uses reduce() not fold().
I would recommend using a backslash instead of ā€œset minusā€ for the solver operator, because ASCII characters are much easier to type.

4 Likes

To add, itā€™s part of CustomStringConvertible protocol.

3 Likes

Hi Nevin, thank you looking at the package! I'll try to address your points:

  • The package is specifically for numerical linear algebra, so I don't think having elements totally unconstrained buys much value. For example, having at least a zero element requirement is useful, for example when dealing with block matrices or for the standard constructor Matrix(rows: Int, columns: Int) . Indeed, I would go as far as constraining matrix elements to be LANumeric , but that would exclude block matrices and integer matrices, which might be nice to have at some point.
  • Similar for AdjointMatrix and TransposeMatrix . In fact, in practical code I would expect people not to use adjoint or transpose much, but using combined operators as ā€²* instead.
  • There are some methods in LANumeric which are supposed to be implemented directly via BLAS, LAPACK and VDSP. It is good for those names to be as close to the original name as possible. It might make sense to separate out these functions (which are the functions implemented in LANumeric+Double , LANumeric+Float , LANumeric+Complex ) into another protocol LANumericCore that is internal to the package. But then again, I think it makes sense to expose them to experienced users who need a broader API surface and want to be generic.
  • The use of toString is quite adhoc and specific to the package. Indeed, description is implemented via toString for matrices in the package. But I cannot change description for complex numbers or vectors as they are external to the package. I also consider toString not as a type conversion but rather as a convenient functionality.
  • It is probably a good idea to rename fold to reduce . I switched forth and back between the names before. I like fold better, reduce has for me somehow the connotation of its parameter function having two arguments of the same type, but that doesn't seem to be so in Swift.
  • I would love to use the backslash instead! It is not possible in Swift though. I agree that ASCII characters are easier to type.

I maintain my viewpoint.

Matrix is a generally useful type for 2D grids of data, regardless of what that data may be.

Constrained extensions make it so that the available functionality matches the capabilities of the Element.

.zero and + come from AdditiveArithmetic.
* comes from Numeric.
etc.

People doing linear algebra would see exactly the same interface with my system as with yours.

I think you are misunderstanding the benefits here.

With separate types, the ā€²* operator could be removed entirely, because uā€² * v would be just as efficient.

Specifically, uā€² would create a new instance that simply wraps a Matrix, with different definitions for its operations.

This is very similar to how someArray.reversed() works.

The names really should follow Swift conventions regardless.

And perhaps they should also be namespaced under .vDSP, .lapack, and .blas, similar to how String has .utf8 and .utf16 views.

Iā€™m just saying, if you want to convert something to a String, then write a String initializer (or a custom interpolation pattern).

Iā€™m not debating the merits of the name, just saying Swift has made a decision.

Ah right, my mistake.

I guess we just have to agree to disagree on most of your points.

But I should be more clear about what I mean by that: In numerical computations, creating an extra memory instance for wrapping is not really a good idea performance-wise. Why wrap something, if you can get rid of the performance penalty in an easy way? But I have to admit that I made a lot of compromises in that regard already in the library, so it wouldn't be that bad, I think it is just not worth the effort.

With regard to naming conventions for LAPACK, BLAS etc. These things have big books of documentation themselves, I don't intend to rewrite them by inventing new names, as swifty as they may be. I would be interested though in encapsulating them somehow so the normal user is ideally not even aware of them.

I don't think it is a good idea to pollute the String namespace with something specific to Matrices. That's just my taste.

As I said, renaming fold to reduce is a good idea.

Itā€™s nice that youā€™re doing this.

But if you donā€™t want to follow Swiftā€™s conventions, I donā€™t really see the point in making a Swift library at all.

5 Likes

It should be very cheap for your Matrix. Even copying an entire Matrix and put it in Adjoint matrix, should cost as much as 1 Array and 2 Int. Copying array is just as expensive as copying a single reference. Itā€™s the mutation thatā€™s expensive (since Array needs to do CoW). Unless weā€™re concerning about arc traffic, which I donā€™t think should be a problem at this scale.

And since the single-value wrapper shares the layout with wrapped object. Maybe (just maybe, I havenā€™t checked), it could avoid copying Matrix altogether.

2 Likes

Yeah, I don't really need comments like that. I follow Swift's conventions where it makes sense, and I don't where it doesn't.

I am using the library and making money by using it. You don't have to.

1 Like

I expect a pull request then soon.

I think you have a fundamental misunderstanding about how memory works in Swift.

Creating the adjoint or transpose would allocate no additional memory.

1 Like

That might be true. Pull request?

Needless to say, itā€™s best to benchmark. But if thereā€™s a regression in doing something similar to this, itā€™s most likely a bug (whether on the programmer side or compiler side).

I think there are two different ways to view matrices. One is as a tool for numerical linear algebra. The current implementation is OK for that purpose.

The other way is to see matrices as a general data structure useful for all sorts of things. That way of viewing matrices is justified, and if somebody wants to extend the library in that way, that might be a good idea. Personally, I do not need it at the moment.

Anyhow, there are tradeoffs between combining operators and creating intermediate types. Since both would result in combinatorial explosions (having myriad * and ' combinations vs having myriad * type overloads) and it affects the ergonomics, so I'm not one to make judgement call. But memory alloc performance shouldn't be the main consideration (they're similar).

I am actually warming up to the idea of having various matrix representations, so that sub matrices etc. would also be represented by views. In the medium term, this is probably a good idea.

1 Like

This is addressed in several places in the Swift Numerics repo, but the taxicab / 1-norm isn't a good option because (like the Euclidean / 2-norm) it can overflow even for finite inputs. The max / āˆž-norm is the only standard norm that's guaranteed to be representable for any finite input. Since all norms on finite-dimensional vector spaces are compatible, this isn't a big deal in general.

I'm actually not totally sure why BLAS chose to use the 1-norm for I*AMAX; the algorithms in question all work just as well with the āˆž-normĀ¹ (and using the āˆž-norm would actually allowed the complex implementation to piggyback on the optimized real implementation when stride == 1), unfortunately, the early BLAS didn't provide a lot of rationale documentation and some of the folks involved are now dead, but I'll ask Velvel and Beresford next time I see them; they're as likely to know as anyone is.

Anyway, long digression, but Swift Numerics will keep the current definition of magnitude, though I could be convinced to also add a binding for the taxicab norm.

Ā¹ the 2-norm would be better algorithmically in some cases, but was historically substantially worse computationally.

2 Likes

It should be possible to implement finding the largest element also for the inf-norm, by running the finding function twice with a stride of 2, so that shouldn't be a problem.

Totally forgot to mention. You have a few of these:

fileprivate func recast<U, V>(_ ptr : UnsafeMutablePointer<U>) -> UnsafeMutablePointer<V> {
    let p = UnsafeMutableRawPointer(ptr)
    return p.assumingMemoryBound(to: V.self)
}

This would be undefined (I think?). You probably want to use one of these: unsafeBitCast(_:to:), numericCast(_:), Int*.init(bitPattern:).

1 Like

I don't think that would be undefined, why? It's the same pointer.