LU decomposing matrix and immutability


I have been trying to create a method for LU decomposing a matrix in Elixir using Doolittle’s algorithm. However, it is defined in an imperative way and trying to implement it functionally has brought about some issues.

The way that feels cleanest and that I’ve worked the most with is to first generate an index set with

#m: height
#n:  width
indexes = for i <- 1..m, do: for j <- 1..n, do: {i, j}

And then map that set to the calculated values, &, fn
  1, j -> a[1][j]          #Guard for the case where k <- 1..0
  i, 1 -> a[i][1]/u[1][1]  #Guard for the case where k <- 1..0
  i, j -> (a[i][j] - Enum.sum(for k <- 1..j-1, do: u[k][j] * l[i][k]))/u[j][j]
  i, j -> a[i][j] - Enum.sum(for k <- 1..i-1, do: u[k][j] * l[i][k])

However, by mapping the values they are all made available at the same time. That is a problem since the values at higher indexes in the matrix depends on the values in the lower ones,

I have thought about implementing it with recursive function calls that updates one element at the time, but then the matrix would need to be copied for every change made. Maybe that is what it takes, but I would love some input! If I have provided to few details, please ask and I will provide more!

Thanks for reading
Sebastian Callh

I would probably put the matrix into a map like enum.into [{{1,1}, 5}, {{1,2}, 7}, {{2,1}, 8}, {{2,2}, 9}], %{} and have 2 matrix = Enum.reduce(indexes, matrix, …) blocks, one for the lower fields and one for the upper fields. Given that you can map the new values into the original A matrix, this should work fine. You could also use a Enum.reduce for the inner ‘p’ loop too.

With the matrix in a Map, you can update the map with Map.put(matrix, {i,j}, value) …

There might be a way to further optimize, but I don’t know the algorithm well enough to suggest anything further.


1 Like

You might want to take a look at tensor, which is a library for sparse matrices in Elixir. it stores matrices internally as maps, and while the Doolittle’s algorithm is not used inside, similar methods are used for recursing over the different fields of the data structures.


As a very minor note, you can do

for i <- 1..m, j <- 1..n, do: {i, j}

instead of nesting for comprehensions.


Hello Steve!

Putting the matrix in a map is something I thought of, for constant access times, so I’m with you on that one. I do not understand how two Enum.reduce(…) would work though, since the first reduce would require values calculated by the second and vice versa, no?



Thanks for the heads-up! I want to implement this myself this time, but I might very well want to use the library in the future!


Hi Ben!

Thanks for that sweet syntactic sugar!


No problem! Do note that it isn’t actually syntactic sugar, but a key part of how for works. It receives a list of generators as well as a list of filters. So for example you can do:

for x <- 1..10, x %2 == 0, y <- 1..x, do: {x, y}

For more information see:

I see! That’s really neat!