Repository: coreylowman/dfdx (language: rust)

Issue:
Tiny loss discrepancy between different matmul features
Might be too low-level/architecture dependant or simply some rounding taking place somewhere, but documenting it here for tracking:

- AMD no features:
```
Epoch 0 in 7.7529317s (241.844 batches/s): avg sample loss 7.75558
Epoch 1 in 7.5624689s (247.935 batches/s): avg sample loss 3.03752
```

- AMD `threaded-cpu`:
```
Epoch 0 in 3.7411362s (501.185 batches/s): avg sample loss 7.76261
Epoch 1 in 3.5128473s (533.755 batches/s): avg sample loss 3.00306
```

- Intel `intel-mkl`:
```
Epoch 0 in 4.0160619s (466.875 batches/s): avg sample loss 7.76278
Epoch 1 in 3.695467s (507.378 batches/s): avg sample loss 3.01809
```

[Relevant matrixmultiply code](https://github.com/bluss/matrixmultiply/blob/master/src/threading.rs)


Discussion / hints:
When manually setting the cache-sensitive MC, KC and NC in the sgemm/dgemm algorithms, I can sort of reproduce the losses:

```
$Env:MATMUL_SGEMM_MC=16;  $Env:MATMUL_SGEMM_KC=16;  $Env:MATMUL_SGEMM_NC=64; cargo run --release --features threaded-cpu --example 06-mnist
Epoch 0 in 5.8709336s (319.370 batches/s): avg sample loss 7.75640
Epoch 1 in 5.7994881s (323.304 batches/s): avg sample loss 3.03374
```

which is like the one with no feature, whereas

```
$Env:MATMUL_SGEMM_MC=128;  $Env:MATMUL_SGEMM_KC=256;  $Env:MATMUL_SGEMM_NC=256; cargo run --release --features threaded-cpu --example 06-mnist
Epoch 0 in 3.7337178s (502.180 batches/s): avg sample loss 7.76261
Epoch 1 in 3.5440562s (529.055 batches/s): avg sample loss 3.00306
```
which is identical to the case with threaded on. So maybe each feature sets the MC, KC, NC variables differently.

Don't know exactly why this has an influence in the result of the matmul though.
Thanks for the additional details/detective work, this is super interesting. It might be worth raising an issue on the matrixmultiply repo just for more info.
Glad that I am not the only one that is curious about/finds it interesting :) Yeah, I plan to open one, just wanted to get some more grasp of the details of blas/blis before I open the issue in matrixmultiply, maybe this would already clarify things. Don't see it having too much time to invest in it (I 'd much rather invest it in the QoL errors issue), but I might wait until the weekend, before I open raise the issue there.
Maybe we should consider enabling std for matrixmultiply by removing `default-features =false` from [here](https://github.com/coreylowman/dfdx/blob/main/Cargo.toml#L31)?

The correct thing to do might be to run some accuracy and speed benchmarks, but bluss suggested that std is recommended. right?
we'll still need it optional for no-std, but we should be able to add it to
```toml
std = [..., "matrixmultiply?/std"]
```
Maybe compare matmul std vs threaded-cpu after that kind of change, and it's possible it will close the performance gap. Then it would indicate matmul's threading doesn't kick in (matrices smaller than our threshold, then). If there's something to gain from threading but we don't do it, maybe our heuristics for how to use threads needs updating.
> Maybe compare matmul std vs threaded-cpu after that kind of change, and it's possible it will close the performance gap. Then it would indicate matmul's threading doesn't kick in (matrices smaller than our threshold, then). If there's something to gain from threading but we don't do it, maybe our heuristics for how to use threads needs updating.

The names of the features have changed since I opened the issue. 

```
[features]
default = ["std", "fast-alloc", "cpu-par-matmul"]
nightly = []

std = ["cudarc?/std"]
fast-alloc = ["std"]
no-std = ["no-std-compat", "dep:spin", "cudarc?/no-std"]

cpu-seq-matmul = ["dep:matrixmultiply"]
cpu-par-matmul = ["std", "dep:matrixmultiply", "matrixmultiply?/threading"]
cpu-mkl-matmul = ["dep:cblas-sys", "dep:libc"]
```

For the tests, I removed from `default` the `cpu-par-matmul` which contains `matrixmultiply?/threading`, so it automatically activates std:

```
cargo run --release --features cpu-seq-matmul --example 06-mnist -- .\examples\tmp\

Epoch 0 in 7.7688808s (241.348 batches/s): avg sample loss 7.75558
Epoch 1 in 7.5353707s (248.827 batches/s): avg sample loss 3.03752
```

```
$Env:MATMUL_NUM_THREADS=1; cargo run --release --features cpu-par-matmul --example 06-mnist -- .\examples\tmp\

Epoch 0 in 4.5212268s (414.710 batches/s): avg sample loss 7.76261
Epoch 1 in 4.640334s (404.066 batches/s): avg sample loss 3.00306
```

```
$Env:MATMUL_NUM_THREADS=4; cargo run --release --features cpu-par-matmul --example 06-mnist -- .\examples\tmp\

Epoch 0 in 3.708347s (505.616 batches/s): avg sample loss 7.76261
Epoch 1 in 3.6519572s (513.423 batches/s): avg sample loss 3.00306
```

Loss discrepancies are like we previously observed them (now explainable). Did you mean ` MATMUL_NUM_THREADS=1`  with `matmul std` though? `cpu-par-matmul` 1 thread is still faster than  `cpu-seq-matmul`, but probably this is just fma vs sse2 because of std. By the way, now the MC, KC, NC constants seem to not influence the result now. I don't know if the combination of features in `threaded-cpu` made it possible back then.
> we'll still need it optional for no-std, but we should be able to add it to
> 
> ```toml
> std = [..., "matrixmultiply?/std"]
> ```

Right now, `matrixmultiply?/std` is always activated from `default`, because `matrixmultiply?/threading` activates it (if I am not misinterpreting something). Any suggestions if you would like to leave it as is or change something inside `default`?
Great, I was just curious about the single threaded performance with std.

> By the way, now the MC, KC, NC constants seem to not influence the result now

They still can, according to my testing. Changing these can have the effect of trading off rounding errors in one place vs somewhere else. (That's not the intended use of these, but a thing they cause.)
> Right now, matrixmultiply?/std is always activated from default, because matrixmultiply?/threading activates it (if I am not misinterpreting something). Any suggestions if you would like to leave it as is or change something inside default?

If `cpu-par-matmul` is activated, then we get `matrixmultiply?/std` due to `matrixmultiply?/threading`.

The issue is if we use `cpu-seq-matmul` (and no par-matmul), we still want `matrixmultiply?/std`. However, we need to support cpu-seq-matmul for no-std, so we can't make cpu-seq-matmul depend on std.
It sounds like adding this will solve that? `std = ["matrixmultiply?/std"]`
Still trying to wrap my head around this. Right now:

```
[features]
default = ["std", "fast-alloc", "cpu-par-matmul"]

std = ["cudarc?/std"]
fast-alloc = ["std"]
no-std = ["no-std-compat", "dep:spin", "cudarc?/no-std"]

cpu-seq-matmul = ["dep:matrixmultiply"]
cpu-par-matmul = ["std", "dep:matrixmultiply", "matrixmultiply?/threading"]
cpu-mkl-matmul = ["dep:cblas-sys", "dep:libc"]
```

To run `no-std` now, we need to run `cargo run --no-default-features --features no-std  <stuff>`, because default has triple `std`, right? This is also forced from the implementation inside `lib.rs`


```
#[cfg(all(feature = "std", feature = "no-std"))]
compile_error!("Can't enable both std and no-std. Set default-features = false to disable std");
``` 

Moreover, `--no-default-features --features no-std`  possibly compiles without problems on CUDA (have no NVIDIA gpu at home, so can't tell), but doesn't compile at all on CPU, everything like `Vec, String, format` etc need `std` when the device is CPU. So since `cpu-seq-matmul` is all about CPU I guess, I am a little confused now.

In any case, this has been a chance for me to learn better the dependencies in rust, especially what the `"package-name?/feature-name"` syntax means.

This also caught my eye in the rust documentation, which is not the way that dfdx handles this, but to be honest no idea how important this is, just making my wall of text larger I guess :) :

> For example, if you want to optionally support [no_std](https://doc.rust-lang.org/reference/names/preludes.html#the-no_std-attribute) environments, do not use a no_std feature. Instead, use a std feature that enables std. For example:
> 
> 
> #![no_std]
> 
> #[cfg(feature = "std")]
> extern crate std;
> 
> #[cfg(feature = "std")]
> pub fn function_that_requires_std() {
>     // ...
> }
  
  
  Could be talking nonsense though, as I said, still trying to wrap my head around this. So ready to implement `std = [..., "matrixmultiply?/std"]` if you are sure it does what you had in mind.
> Moreover, --no-default-features --features no-std possibly compiles without problems on CUDA (have no NVIDIA gpu at home, so can't tell), but doesn't compile at all on CPU, everything like Vec, String, format etc need std when the device is CPU. So since cpu-seq-matmul is all about CPU I guess, I am a little confused now.

Ahh good catch, looks like some recent PRs broke this. Looks like we'll need to add a workflow for monitoring no-std support.

> This also caught my eye in the rust documentation, which is not the way that dfdx handles this, but to be honest no idea how important this is, just making my wall of text larger I guess :) :

Yeah I recently moved away from this. Basically no-std was adding a couple dependencies (no-std-compat, and spin) to ALL std builds. This increases compilation time just to support the no-std case. So even folks that don't care about no_std would pay for compilation times.

> So ready to implement std = [..., "matrixmultiply?/std"] if you are sure it does what you had in mind.

Yep lets go for it!

If you want to be sure, you `cargo tree` is a nifty tool. You can run

- threaded matrixmultiply case `cargo tree`
- non threaded matrixmultiply case `cargo tree --no-default-features -F std,cpu-seq-matmul`
- no matrixmultiply case `cargo tree --no-default-features -F std`
- - no std case `cargo tree --no-default-features -F no-std`

Thanks for the `cargo tree` tip, one more thing I get the chance to learn! Will test it tomorrow (not at the computer atm).

My main point in all the above was that there might be no any practical difference in the change, as no `cpu-matmul` looked relevant to me in a `no-std` scenario, so we might as well leave it as is.

If you guide me a bit what you would like to test in the no-std workflow, I could work on that as well (only condition would be that no cuda is necessary for its development, as I have no nvidia gpu at home).

If you prefer to take care of the workflow yourself later at some point, I will pack the `std = [..., "matrixmultiply?/std"]` change in the PR with the QoL Errors for mismatched dims, so that you don't have to merge 2 times 👍  


> If you guide me a bit what you would like to test in the no-std workflow, I could work on that as well (only condition would be that no cuda is necessary for its development, as I have no nvidia gpu at home).

Yeah no cuda necessary! We have someone fixing the no-std stuff now (#615), and adding a workflow to check it automatically when you open a PR.

After that's merged in we can just run `cargo test --tests` in all the scenarios I laid out above

> If you prefer to take care of the workflow yourself later at some point, I will pack the std = [..., "matrixmultiply?/std"] change in the PR with the QoL Errors for mismatched dims, so that you don't have to merge 2 times 👍

Let's keep them separate, merging is cheap and it's good to separate distinct fixes/features 😄 

Fix patch:
diff --git a/Cargo.toml b/Cargo.toml
index 26920ee6a..61a807a44 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -49,7 +49,7 @@ glob = { version = "0.3.1", optional = true }
 default = ["std", "fast-alloc", "cpu-par-matmul"]
 nightly = []
 
-std = ["cudarc?/std"]
+std = ["cudarc?/std", "matrixmultiply?/std"]
 fast-alloc = ["std"]
 no-std = ["no-std-compat", "dep:spin", "cudarc?/no-std"]
 


Test patch: