I tried to post a question on https://groups.google.com/forum/#!forum/uncomplicate but it did not seem to go through. So I will ask it here. I have been following along with the linear algebra refresher posts - but also tried to use neantherthal to solve the systems of linear equations (using sv!). And the answers were different from those in the book. I found that by changing the order of the rows I could get the right answer. I first thought this was because of the underlying lapack implementation but when I used python numpy I always got the same answer as the book independent of the order of the rows.

Here are two examples:

(let [cs (dge 3 1 [16 8 0])]
  (sv! (dge 3 3 [0 4 1
                 4 0 1
                 1 1 -1])
       cs)
  cs)

gives:

#RealGEMatrix[double, mxn:3x1, order:column, offset:0, ld:3]

  1.00 
  3.00 
  4.00 

but:

(let [cs (dge 3 1 [0 8 16])]
  (sv! (dge 3 3 [1 1 -1
                 4 0 1
                 0 4 1])
       cs)
  cs)

gives:

#RealGEMatrix[double, mxn:3x1, order:column, offset:0, ld:3]

 -9.33 
  2.33 
  4.33 

I guess my question is: “Why am I getting different answers?”

@rrottier You swapped the variables instead of the equations in your second example, because you forgot that the default matrix layout is column-oriented. Here's what you intended to do (the answer is also what you'd expect):

(let [cs (dge 3 1 [0 16 8])]
   (sv! (dge 3 3 [1 0 4 
                          1 4 0 
                          -1 1 1])
      cs)
   cs)

the answer is:

#RealGEMatrix[double, mxn:3x1, order:column, offset:0, ld:3]

   1.00 
   3.00 
   4.00 

Now, maybe it is inconvenient to you to work with this "rotated" perspective. No problem, just add the option {:order :row} to dge and the matrix will be row-oriented:

(let [cs (dge 3 1 [0 8 16] {:order :row})]
  (sv! (dge 3 3 [1 1 -1
                 4 0 1
                 0 4 1]
            {:order :row})
       cs)
  cs)

=> #RealGEMatrix[double, mxn:3x1, order:row, offset:0, ld:1]

   1.00 
   3.00 
   4.00 

@blueberry Thanks. I did not realise the default matrix layout is column oriented. Thanks for pointing that out.

Is this a property of the underlying lapack or is it a neantherthal design decision?

@rrottier It's like the parentheses - after a few hours, you will not notice it :) BTW. matrix string prints out that information, and if I remember well, I mention it quite a few times in the tutorials. Please do this kind of feedback, and, if possible, suggest how to make it more obvious.

both

PS. Do I really need to specify {:order :row} on (dge 3 1 ...)? It seems like as it is only one column the orientation should be obvious.

Yeah #RealGEMatrix[double, mxn:3x1, order:column, offset:0, ld:3] I did not really parse the type signature, was too focussed on the values 😊

Both LAPACK and neanderthal handle both, but column is the default in literature. In lapack, column is the only option for some operations, but neanderthal finds a way to provide roe even in most of these cases.

“column is the default in literature” I thought that row is the default in literature.

You need it even for 3x1 because the LAPACK solver requires matching layout for both matrices.

At least the literature that i read for systems of linear equations always represent the 2x + 3y + 4z as [2 3 4]

but it is one vector :) that does not have order

order is a computing detail. your books are math textbooks.

Ah… that is very true - I am only now stumbling into the computing stuff.

the first rule of numerical computing is: math book is not enough :)

I thought functional programming’s promise is: “You can write it just like the math”

if you ignore that, you'll be heavily bitten by floating point arythmetic.

I do not know about that promise, but It is not possible to fullfil

You can write it whatever you want, provided that you do not need to compute it :)

Numpy seems to do some transparent casting then…

Because it behaves as row oriented.

sorry I should have said “non-transparent”

Neanderthal also supports rows and columns. It's just about the default order, and does not have anything to do with floating point precision. In computing, a*b is not always equal to b*a.

I think I just stumbled into the column orientation with the mm function as well. I tried to mutiply: (dge 2 2 [1 2 3 4]) with (dge 2 1 [1 2]) and got a dimension error. Does that make sense?

I could also implement such castings, but that is something i want to avoid since it can shoot you in the foot when you need performance.

no

Sorry - just edited it. Got it the wrong way round

that does not ave anything to do with order. 2x1 * 2x2 is not allowed by math rules

btw. neanderthal supports mixing orders in most operations. just not in solvers.

RE: The castings - the way you do it is probably better, I agree that you want to be transparent to the underlying implementation. It is just confusing coming form another abstraction

I was trying to do this in neantherthal.

post the code an i'll tell you where you made a mistake

this works

(let [a (dge 2 2 [1 2 3 4])
      b (dge 1 2 [2 5])]
  (mm a b))

you are trying to multiply 2x2 and 1x2

should be 2 1, not 2 w

not 1 2

Sorry same mistake as previous

(let [a (dge 2 2 [1 2 3 4])
      b (dge 2 1 [2 5])]
  (mm a b))

column orientation does not require you to rotate the whole world :)

Now it works but gives the wrong answer

[17
 24]

but if the matrix is column oriented, it transfers the 1d vector [1 2 3 4] into columns.

it should be [1 3 2 4]

My head hurts…

look, the only thing you need is this:

Ok, I think I understand where all my issues are coming from… I just need to change the way I think about it.

matrix is 2d, vector is 1d

you wand to create a matrix by supplying the data in a vector

read my latest blog post. i wrote the whole blog post about these issues

I actually read it.

Did not understand the implications of the column orientation though.

you probably do not need to change the way you think about it, you just need to be aware of a few etails.

the use matrices efficiently one?

But I get the reasons - in haskell I started playing with dense which uses the same kind of abstraction.

http://dragan.rocks/articles/17/Clojure-Numerics-1-Use-Matrices-Efficiently

Thanks for the help.

What I want, is to lay AA in memory in a way that is efficient, but still easy to work with. Obviously, those two dimensions have to be projected into one. But how? Should it be [1112|2313|1−1−2−6][1112|2313|1−1−2−6], or [121|13−1|11−2|23−6][121|13−1|11−2|23−6]? The answer is: it depends on how your algorithm accesses it most often. Neanderthal gives you both options. When you create any kind of matrix, you can specify whether you want it to be column-oriented (:column, which is the default), or row-oriented (:row). In the following example, we will use CPU matrices from the native namespace. The same options also work for functions that create GPU CUDA matrices (cuda namespace), or OpenCL's GPU and CPU matrices (opencl namespace).

and after that...

Clojure sequence that we used as data source has been read column-by-column.

“which is the default” is what I missed.

🙂

But this conversation gives me material for a nice blog post, so thank you.

“No worries” as they say in Australia. Thanks for the linear algebra series.

I really enjoy it.