Rank one (outer) product

[ghbrown]

Description

A simple function to compute the rank one (a.k.a. outer) product of two vectors, which would return a matrix.

I believe the current way to do this (or at least the way I do) involves explicit looping or reshapeing and then matmul (even worse than looping probably). I just suggest the mathematics be hidden away from a user inside the function, so they can pass two true vectors in and it will return a matrix.

Rank one updates are sufficiently common that I think we should include a specific function to accommodate this. Here is my crack at a simple "column-major-optimized" version (I'd do speed tests comparing the various methods if we want this function):

function rop(u,v) result(mat)
    implicit none
    real, intent(in) :: u(:), v(:)
    real :: mat(size(u),size(v))
    integer :: col
    do col=1,size(v)
        mat(:,col) = v(col) * u
    end do
end function

Options/alternatives to a dedicated function:

• this functionality may eventually be provided by a more general tensor product routine (like those below), so perhaps this would just be a holdover
• something like this likely exists in BLAS somewhere in which case this function would be a wrapper for that if desired
• there could also be a good reason (that I don't) know why this feature doesn't exist, hasn't been proposed, and shouldn't be implemented (perhaps like the fact it requires quadratic storage for a linear amount of data in which case it might be preferable to update all elements in a loop instead, who knows?)

Prior Art

Numpy does not provide a rank-one-product-specific function, but provides this functionality through multiple functions, including:

• numpy.outer(u,v)
• numpy.tensordot(u,v,axes=0)
• numpy.einsum('i,j',u,v) or numpy.einsum('i,j->ij',a,a)
  given vectors u, v.

Relevant but separate discussion
Perhaps this is just a result of my field of interest, but I'd like to bring up the discussion of terminology for the number of array dimensions (only since it's relevant here):

• I believe the current Fortran terminology is "rank", this is not optimal in my eyes because rank has a conflicting linear algebra definition (especially problematic since Fortran is commonly used for linear algebra); Python/Numpy also uses rank, but that doesn't make it "right"
  -"order", "mode", or "way" are my suggestions and what many people use as an alternative to rank in literature on multi-way arrays beyond matrices
• "dimension" seems to make sense but is almost ubiquitously used to describe the length of a mode, so I would say it's out of the question

[ivan-pi]

Here is an implementation I have used previously:

module outer_prod_mod
  implicit none
  integer, parameter :: wp = kind(1.0d0)
contains
!>
!!   Outer product of two vectors
!!
  pure function outer_prod(u,v) result(p)
    real(wp), intent(in) :: u(:), v(:)
    real(wp) :: p(size(u),size(v))

    integer :: n, m

    n = size(u)
    m = size(v)
    p = spread(u,dim=2,ncopies=m)*spread(v,dim=1,ncopies=n)
  end function
end module

[ivan-pi]

We also have an issue open for einsum (#376). But I'd be happy to have both just like NumPy has. The einsum notation is more suitable for complex tensor reductions.

Maybe we could name it outer_product, since Fortran already has dot_product. If you want a shorter name you can always abbreviate at import, e.g.: use stdlib_linalg, only: outer => outer_product.

[ghbrown]

I agree with having both for convenience, and the outer_product name is good with me.

I'm happy to work on this as an intro to the stdlib project if you'd like. I'm all set up on the software side and have read the workflow, style guide, and CONTRIBUTING.md. However, I'm not familiar with fypp, and haven't contributed to an open source project before, so things may be a bit slow at first. Is there perhaps a more detailed guide to the project structure and tools on the fortran-lang Discourse somewhere that could help me get my bearings?

Also, what is the preference (if any) for spread() over other methods for the rank one product? I hadn't hear of it before.

[ejovo13]

Hey GHBrown!

I'm in the same boat as you and just this last week I proposed my first pull request ever on github ever to this project (which im still cleaning up).

Anyways, concerning fypp, you'll pick it up really quick. It's a preprocessor that's called on .fypp files to provide meta-programming tools like generic programming. This is awesome because Fortran has generic procedures but nothing like C++ templates where a new function signature will get created at compile time. So, using fypp python loops that iterate through all the fortran types, we can emulate template programming by creating function and signature definitions for all reasonable combinations. This extends into creating generic functions on arrays with varying ranks, say up until 15-rank inputs.

My advice for learning how to use fypp is to read the docs, start with the introduction and definitely check out the generic programming paragraph that I linked you. Afterwards, see how fypp is implemented in the stdlib source folder. Be sure to read through common.fypp to see some common variables that are used throughout the source code, like REAL_KINDS_TYPES. If you have experience in python, you'll be able to use it really easily. If not, you could try what I did, which was to open up a jupyter notebook and just play around with the variables are defined in common.fypp and how you can access elements of the lists/tuples, especially when looping.

Basically what happens is during the build process the source code files are treated by fypp (at least the ones that you designate in the cmake file) and then the compiler will run through the processed files. You can see the outputs of fypp in your build/src folder when building with cmake.

The best teacher is observation - I've noticed that the project structure for source code is defining all of your interfaces in a top level module (like stdlib_math for example) and then placing your function's source code in a submodule in the file stdlib_math_subm.fypp. Just poke around the source code a little bit, and imitate how other authors have implemented fypp.

Once you've created your source files, make sure you add them to the src/CMakeLists.txt file, and also to the Makefile.manual file.

Good luck with your first contribution to open source and let me know if you have other questions

[ivan-pi]

Also, what is the preference (if any) for spread() over other methods for the rank one product? I hadn't hear of it before.

The only preference was the implementation is very short. Leaving out the keywords the main body of the routine would be:

n = size(u)
m = size(v)
p = spread(u,2,m)*spread(v,1,n)

The function spread will replicate a source array ncopies times along a specified dimension dim. It is similar to numpy.tile. Upon more thought it might not be the most efficient way since it is likely the compiler will create two temporary arrays, which could become expensive if n and m are large. I think it's okay if we start with your implementation.

@ejovo13 has given a great introduction to fypp. For the beginning you can just start with the default real precision. Once you are satisfied, you will just need to add a preprocessor loop for other types. It is best if this is done in a submodule to keep files small improving both compilation time and memory usage.

To prepare your contribution you should first create a fork of stdlib and clone it to your local PC. Afterwards make a new branch (e.g. git checkout -b outer_product). Then you can add the code necessary to stdlib_math and also some tests in a new file test_outer_product.f90. At that point you can probably already create a draft pull request to solicit help and get some reviews.

[ghbrown]

Thank you both for the advice and information! That cleared up a lot of the big picture questions I had, and I'll consult the docs for the rest.

@ivan-pi, would you like the interface in stdlib_math or stdlib_linalg? You mentioned math, but I was expecting you to suggest linalg.

[ivan-pi]

...but I was expecting you to suggest linalg.

You are right, linalg is more suitable.

[ghbrown]

I'm going to close this issue. If this is premature, please feel free to reopen.