# What is so different in each of these matrix operations?

**URL:** <https://discourse.julialang.org/t/what-is-so-different-in-each-of-these-matrix-operations/78637>\
**Category:** General Usage\
**Tags:** question\
**Created:** [March 28, 2022, 9:13pm UTC](https://discourse.julialang.org/t/what-is-so-different-in-each-of-these-matrix-operations/78637 "2022-03-28T21:13:05Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![csubhodeep](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/csubhodeep/32/34967_2.png) [@csubhodeep](https://discourse.julialang.org/u/csubhodeep)\
**Post date:** [March 28, 2022, 9:13pm UTC](https://discourse.julialang.org/t/what-is-so-different-in-each-of-these-matrix-operations/78637/1 "2022-03-28T21:13:05Z")

</div>

Hello all,

I am fairly new to Julia and have just started learning the basics of `Matrix`. I wanted to get a brief understanding of what is so different under the hood that is causing the huge time differences in each of the following approaches when technically they are basically the same operations yielding the same result.

The functions:

```julia
using LinearAlgebra

function train_func(x_train::Matrix{Float64}, y_train::Matrix{Float64})::Matrix{Float64}

    return pinv(x_train)*y_train

end;

function train_func_2(x_train::Matrix{Float64}, y_train::Matrix{Float64})::Matrix{Float64}

    # this is assuming n_rows > n_columns
    return inv(transpose(x_train)*x_train)*transpose(x_train)*y_train
    
end;

function train_func_3(x_train::Matrix{Float64}, y_train::Matrix{Float64})::Matrix{Float64}

    return x_train\y_train

end;

```

Testing with a toy example

```julia
julia> x = [1. 2.; 3. 4.; 5. 6.]
3×2 Matrix{Float64}:
 1.0 2.0
 3.0 4.0
 5.0 6.0

julia> y = reshape([1. 2. 3.], 3, 1)
3×1 Matrix{Float64}:
 1.0
 2.0
 3.0

julia> train_func(x,y) ≈ train_func_2(x,y) ≈ train_func_3(x,y)
true

```

Now if we do the same with slightly larger matrices

```julia
julia> n, m = 20000, 1000
(20000, 1000)

julia> x = rand(Int64).*rand(Float64, (n,m));

julia> y = rand(Float64, (n,1));

julia> @time train_func(x,y); @time train_func_2(x,y); @time train_func_3(x,y);
  3.946035 seconds (30 allocations: 648.657 MiB, 0.31% gc time)
  0.247796 seconds (9 allocations: 15.771 MiB)
  7.486294 seconds (6.03 k allocations: 168.779 MiB, 0.51% gc time)

```

Runtime details:

```julia
julia> versioninfo()
Julia Version 1.7.2
Commit bf53498635 (2022-02-06 15:21 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i5-10210U CPU @ 1.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, skylake)

```

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [March 28, 2022, 9:51pm UTC](https://discourse.julialang.org/t/what-is-so-different-in-each-of-these-matrix-operations/78637/2 "2022-03-28T21:51:11Z")

</div>

Oops, that looks like a terrible counter example to what we usually trust…

```julia
using LinearAlgebra
using BenchmarkTools

function train_func(x_train, y_train)
    return pinv(x_train)*y_train
end;

function train_func_2(x_train, y_train)
    # this is assuming n_rows > n_columns
    return inv(transpose(x_train)*x_train)*transpose(x_train)*y_train
end;

function train_func_3(x_train, y_train)
    return x_train\y_train
end;

n, m = 20000, 1000
x = rand(Int64).*rand(Float64, (n,m));
y = rand(Float64, n);
result1 = train_func(x,y); 
result2 = train_func_2(x,y); 
@assert result2 ≈ result1
result3 = train_func_3(x,y);
@assert result3 ≈ result2

@btime train_func($x,$y); 
@btime train_func_2($x,$y); 
@btime train_func_3($x,$y);

```

yielding

```julia
  2.713 s (30 allocations: 648.66 MiB)
  194.618 ms (9 allocations: 15.77 MiB)
  4.097 s (6031 allocations: 168.78 MiB)

```

Kudos for bringing this up!

Edit: but luckily this

```julia
function train_func_4(x_train, y_train)
    return (transpose(x_train)*x_train)\(transpose(x_train)*y_train)
end;

```

yielding

```julia
  125.174 ms (7 allocations: 15.28 MiB)

```

repairs my trust again…

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [March 29, 2022, 12:30am UTC](https://discourse.julialang.org/t/what-is-so-different-in-each-of-these-matrix-operations/78637/3 "2022-03-29T00:30:53Z")

</div>

Here you go [Efficient way of doing linear regression - #33 by stevengj](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/33)

---

<div class="post-metadata">

**Author:** ![csubhodeep](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/csubhodeep/32/34967_2.png) [@csubhodeep](https://discourse.julialang.org/u/csubhodeep)\
**Post date:** [March 29, 2022, 7:06am UTC](https://discourse.julialang.org/t/what-is-so-different-in-each-of-these-matrix-operations/78637/4 "2022-03-29T07:06:14Z")

</div>

Thanks a lot @zdenek_hurak for sharing the better implementation and the link to the right thread.