I ported fixest core algorithms to Armadillo (C++)

[pachadotdev]

Hi @s3alfisc

After finding this, I ported the core algorithms from fixest to Armadillo, so these can be used either in Python, R, etc.

I put my version here: https://github.com/pachadotdev/capybara/tree/reghdfe-ideas

While it is an R package, all the link functions, deviance, etc have been written on C++ side, and the R functions only validate the data and pass it to Armadillo for the computation.

I still need to fix a few issues

[apoorvalal]

Hi Pacha, this looks like a cool package. I've looked into porting over Laurent's optimized cpp demeaner in the past but the extensive dependencies on Rcpp (as opposed to vanilla cpp) got in the way; based on a quick glance at src/ in Capybara this shouldn't be a problem?

If so, i think it would be pretty straightforward to provide python bindings via pybind11 or nanobind; I use pybind11 in pyensmallen. Any chance you could take a stab at porting over some of this on a branch to this repo and we can test out having capybara as a submodule we pull from your repo?

[pachadotdev]

hi @apoorvalal !

exactly!

what I did follows an "epoxy glue" approach in two parts:

1. Armadillo wrappers (no Rcpp): https://github.com/pachadotdev/cpp11armadillo
2. Armadillo demeaning, fitting, etc (based on fixest, alpaca, ppmlhdfe, reghdfe and base R): https://github.com/pachadotdev/capybara/blob/reghdfe-ideas/src/03_demean_variables.h

the advantage of this is that what is done on R side is minimal, just some data.table operations

I looked at base R internals to write all the link functions, deviance, etc using Armadillo, so as you mention, plugging this into Python with pybind11 and pandas should be kinda direct

As a side thing, Kendall's correlation for PPML pseudo-R2 is too slow in base R, so I wrote this https://github.com/pachadotdev/kendallknight

You can find more details here:

Vargas Sepúlveda, Mauricio. 2025. “Kendallknight: An R Package for Efficient Implementation of Kendall's Correlation Coefficient Computation.” PLOS ONE 20 (6): e0326090. https://doi.org/10.1371/journal.pone.0326090.

Vargas Sepúlveda, Mauricio and Schneider Malamud, Jonathan. 2025. “cpp11armadillo: An R package to use the Armadillo C++ library.” SoftwareX 30 (May): 102087. https://doi.org/10.1016/j.softx.2025.102087.

[pachadotdev]

ps. let me fix some errors in the new branch that I do not appear with alternating projections in the main branch, and then I'll try to plug this into Python

the errors are:

> test()
ℹ Testing capybara
✔ | F W  S  OK | Context
✖ | 1        7 | apes-bias                                                                              
────────────────────────────────────────────────────────────────────────────────────────────────────────
Failure (test-apes-bias.R:44:3): apes/bias works
round(coef(apes1), 2) (`actual`) not equal to `apes2` (`expected`).

  `actual`: 0.03
`expected`: 0.05
────────────────────────────────────────────────────────────────────────────────────────────────────────
✔ |          5 | autoplot                                                                               
✔ |          1 | deterministic                                                                          
✔ |         25 | errors                                                                                 
✔ |          3 | feglm-atypical-columns                                                                 
✖ | 6        8 | feglm                                                                                  
────────────────────────────────────────────────────────────────────────────────────────────────────────
Failure (test-feglm.R:33:3): feglm is similar to glm
unname(coef(mod_binom) - coef(mod_binom_base)[2:3]) (`actual`) not equal to c(0, 0) (`expected`).

  `actual`: -0.66 0.02
`expected`:  0.00 0.00

Failure (test-feglm.R:108:3): predicted values increase the error outside the inter-quartile range for GLMs
`pred1_binom` (`actual`) not equal to unname(pred1_base_binom) (`expected`).

  `actual`: 0.97 0.79 0.97 0.03 0.07 0.01 0.01 0.02 0.02 0.00 and 7 more...
`expected`: 0.98 0.88 0.95 0.03 0.14 0.01 0.01 0.06 0.06 0.01           ...

Failure (test-feglm.R:109:3): predicted values increase the error outside the inter-quartile range for GLMs
`pred2_binom` (`actual`) not equal to unname(pred2_base_binom) (`expected`).

  `actual[1:6]`: 0.02 0.01 0.02 0.00 0.00 0.00
`expected[1:6]`: 0.06 0.01 0.10 0.00 0.00 0.00

  `actual[7:15]`: 0.99 1.00 1.00 0.18 0.00 1.00 0.99 1.00 0.06
`expected[7:15]`: 0.99 1.00 1.00 0.44 0.01 1.00 0.98 1.00 0.23

Failure (test-feglm.R:139:3): predicted values increase the error outside the inter-quartile range for LMs
`pred1_binom` (`actual`) not equal to unname(pred1_base_binom) (`expected`).

  `actual`: 0.97 0.79 0.97 0.03 0.07 0.01 0.01 0.02 0.02 0.00 and 7 more...
`expected`: 0.98 0.88 0.95 0.03 0.14 0.01 0.01 0.06 0.06 0.01           ...

Failure (test-feglm.R:140:3): predicted values increase the error outside the inter-quartile range for LMs
`pred2_binom` (`actual`) not equal to unname(pred2_base_binom) (`expected`).

  `actual[1:6]`: 0.02 0.01 0.02 0.00 0.00 0.00
`expected[1:6]`: 0.06 0.01 0.10 0.00 0.00 0.00

  `actual[7:15]`: 0.99 1.00 1.00 0.18 0.00 1.00 0.99 1.00 0.06
`expected[7:15]`: 0.99 1.00 1.00 0.44 0.01 1.00 0.98 1.00 0.23

Failure (test-feglm.R:157:3): proportional regressors return NA coefficients
predict(fit2) (`actual`) not equal to unname(predict(fit1)) (`expected`).

  `actual`: -0.07 -0.07 -0.10 -0.09 -0.10 -0.10 -0.07 -0.10 -0.10 -0.06 and 990 more...
`expected`: -0.07 -0.06 -0.08 -0.08 -0.08 -0.08 -0.06 -0.09 -0.09 -0.06             ...
────────────────────────────────────────────────────────────────────────────────────────────────────────
✔ |         27 | felm                                                                                   
✖ | 1        0 | fenegbin                                                                               
────────────────────────────────────────────────────────────────────────────────────────────────────────
Failure (test-fenegbin.R:27:3): fenegbin is similar to fixest
coef(mod_base)[2] (`actual`) not equal to coef(mod)[1] (`expected`).

  `actual`: -0.18
`expected`: -0.19
────────────────────────────────────────────────────────────────────────────────────────────────────────
✔ |          2 | fepoisson-extra                                                                        
✖ | 2       18 | fepoisson                                                                              
────────────────────────────────────────────────────────────────────────────────────────────────────────
Failure (test-fepoisson.R:43:3): fepoisson is similar to fixest
0.05 is not strictly more than `dist_variation`. Difference: -0.0022

Failure (test-fepoisson.R:103:3): fepoisson is similar to fixest
`dist_variation` is not strictly less than 0.05. Difference: 0.0668
────────────────────────────────────────────────────────────────────────────────────────────────────────
✖ | 7        3 | fixed-effects                                                                          
────────────────────────────────────────────────────────────────────────────────────────────────────────
Failure (test-fixed-effects.R:27:3): fixed_effects is similar to glm
rep(0, 10) (`actual`) not equal to c1 - c2 (`expected`).

  `actual`:  0.00  0.00 0.00  0.00 0.00 0.00  0.00  0.00  0.00  0.00
`expected`: -0.01 -0.01 0.01 -0.03 0.02 0.01 -0.02 -0.02 -0.03 -0.01

Failure (test-fixed-effects.R:53:3): fixed_effects is similar to glm
unname(coef(mod_binom) - coef(mod_binom_base)[2:3]) (`actual`) not equal to c(0, 0) (`expected`).

  `actual`: -0.66 0.02
`expected`:  0.00 0.00

Failure (test-fixed-effects.R:54:3): fixed_effects is similar to glm
unname(coef(mod_binom) - coef(mod_binom_fixest)) (`actual`) not equal to c(0, 0) (`expected`).

  `actual`: -0.66 0.02
`expected`:  0.00 0.00

Failure (test-fixed-effects.R:68:3): fixed_effects is similar to glm
fe_binom - fe_binom_base (`actual`) not equal to c(0, 0, 0) (`expected`).

  `actual`: 1.6 1.0 0.8
`expected`: 0.0 0.0 0.0

Failure (test-fixed-effects.R:69:3): fixed_effects is similar to glm
unname(fe_binom - fe_binom_fixest$cyl) (`actual`) not equal to c(0, 0, 0) (`expected`).

  `actual`: 1.6 1.0 0.8
`expected`: 0.0 0.0 0.0

Failure (test-fixed-effects.R:82:3): fixed_effects is similar to glm
fe_gamma - fe_gamma_base_corrected (`actual`) not equal to c(0, 0, 0) (`expected`).

`names(actual)` is a character vector ('(Intercept)', '(Intercept)', '(Intercept)')
`names(expected)` is absent

  `actual`: 0.04 0.04 0.03
`expected`: 0.00 0.00 0.00

Failure (test-fixed-effects.R:95:3): fixed_effects is similar to glm
fe_invgauss - fe_invgauss_base_corrected (`actual`) not equal to c(0, 0, 0) (`expected`).

`names(actual)` is a character vector ('(Intercept)', '(Intercept)', '(Intercept)')
`names(expected)` is absent
────────────────────────────────────────────────────────────────────────────────────────────────────────
✔ |          4 | group-sums                                                                             
✔ |          2 | no-fixed-effects                                                                       
✔ |          3 | vcov                                                                                   
✔ |          6 | weights                                                                                

══ Results ═════════════════════════════════════════════════════════════════════════════════════════════
Duration: 1.3 s

── Failed tests ────────────────────────────────────────────────────────────────────────────────────────
Failure (test-apes-bias.R:44:3): apes/bias works
round(coef(apes1), 2) (`actual`) not equal to `apes2` (`expected`).

  `actual`: 0.03
`expected`: 0.05

Failure (test-feglm.R:33:3): feglm is similar to glm
unname(coef(mod_binom) - coef(mod_binom_base)[2:3]) (`actual`) not equal to c(0, 0) (`expected`).

  `actual`: -0.66 0.02
`expected`:  0.00 0.00

Failure (test-feglm.R:108:3): predicted values increase the error outside the inter-quartile range for GLMs
`pred1_binom` (`actual`) not equal to unname(pred1_base_binom) (`expected`).

  `actual`: 0.97 0.79 0.97 0.03 0.07 0.01 0.01 0.02 0.02 0.00 and 7 more...
`expected`: 0.98 0.88 0.95 0.03 0.14 0.01 0.01 0.06 0.06 0.01           ...

Failure (test-feglm.R:109:3): predicted values increase the error outside the inter-quartile range for GLMs
`pred2_binom` (`actual`) not equal to unname(pred2_base_binom) (`expected`).

  `actual[1:6]`: 0.02 0.01 0.02 0.00 0.00 0.00
`expected[1:6]`: 0.06 0.01 0.10 0.00 0.00 0.00

  `actual[7:15]`: 0.99 1.00 1.00 0.18 0.00 1.00 0.99 1.00 0.06
`expected[7:15]`: 0.99 1.00 1.00 0.44 0.01 1.00 0.98 1.00 0.23

Failure (test-feglm.R:139:3): predicted values increase the error outside the inter-quartile range for LMs
`pred1_binom` (`actual`) not equal to unname(pred1_base_binom) (`expected`).

  `actual`: 0.97 0.79 0.97 0.03 0.07 0.01 0.01 0.02 0.02 0.00 and 7 more...
`expected`: 0.98 0.88 0.95 0.03 0.14 0.01 0.01 0.06 0.06 0.01           ...

Failure (test-feglm.R:140:3): predicted values increase the error outside the inter-quartile range for LMs
`pred2_binom` (`actual`) not equal to unname(pred2_base_binom) (`expected`).

  `actual[1:6]`: 0.02 0.01 0.02 0.00 0.00 0.00
`expected[1:6]`: 0.06 0.01 0.10 0.00 0.00 0.00

  `actual[7:15]`: 0.99 1.00 1.00 0.18 0.00 1.00 0.99 1.00 0.06
`expected[7:15]`: 0.99 1.00 1.00 0.44 0.01 1.00 0.98 1.00 0.23

Failure (test-feglm.R:157:3): proportional regressors return NA coefficients
predict(fit2) (`actual`) not equal to unname(predict(fit1)) (`expected`).

  `actual`: -0.07 -0.07 -0.10 -0.09 -0.10 -0.10 -0.07 -0.10 -0.10 -0.06 and 990 more...
`expected`: -0.07 -0.06 -0.08 -0.08 -0.08 -0.08 -0.06 -0.09 -0.09 -0.06             ...

Failure (test-fenegbin.R:27:3): fenegbin is similar to fixest
coef(mod_base)[2] (`actual`) not equal to coef(mod)[1] (`expected`).

  `actual`: -0.18
`expected`: -0.19

Failure (test-fepoisson.R:43:3): fepoisson is similar to fixest
0.05 is not strictly more than `dist_variation`. Difference: -0.0022

Failure (test-fepoisson.R:103:3): fepoisson is similar to fixest
`dist_variation` is not strictly less than 0.05. Difference: 0.0668

Failure (test-fixed-effects.R:27:3): fixed_effects is similar to glm
rep(0, 10) (`actual`) not equal to c1 - c2 (`expected`).

  `actual`:  0.00  0.00 0.00  0.00 0.00 0.00  0.00  0.00  0.00  0.00
`expected`: -0.01 -0.01 0.01 -0.03 0.02 0.01 -0.02 -0.02 -0.03 -0.01

Failure (test-fixed-effects.R:53:3): fixed_effects is similar to glm
unname(coef(mod_binom) - coef(mod_binom_base)[2:3]) (`actual`) not equal to c(0, 0) (`expected`).

  `actual`: -0.66 0.02
`expected`:  0.00 0.00

Failure (test-fixed-effects.R:54:3): fixed_effects is similar to glm
unname(coef(mod_binom) - coef(mod_binom_fixest)) (`actual`) not equal to c(0, 0) (`expected`).

  `actual`: -0.66 0.02
`expected`:  0.00 0.00

Failure (test-fixed-effects.R:68:3): fixed_effects is similar to glm
fe_binom - fe_binom_base (`actual`) not equal to c(0, 0, 0) (`expected`).

  `actual`: 1.6 1.0 0.8
`expected`: 0.0 0.0 0.0

Failure (test-fixed-effects.R:69:3): fixed_effects is similar to glm
unname(fe_binom - fe_binom_fixest$cyl) (`actual`) not equal to c(0, 0, 0) (`expected`).

  `actual`: 1.6 1.0 0.8
`expected`: 0.0 0.0 0.0

Failure (test-fixed-effects.R:82:3): fixed_effects is similar to glm
fe_gamma - fe_gamma_base_corrected (`actual`) not equal to c(0, 0, 0) (`expected`).

`names(actual)` is a character vector ('(Intercept)', '(Intercept)', '(Intercept)')
`names(expected)` is absent

  `actual`: 0.04 0.04 0.03
`expected`: 0.00 0.00 0.00

Failure (test-fixed-effects.R:95:3): fixed_effects is similar to glm
fe_invgauss - fe_invgauss_base_corrected (`actual`) not equal to c(0, 0, 0) (`expected`).

`names(actual)` is a character vector ('(Intercept)', '(Intercept)', '(Intercept)')
`names(expected)` is absent

[ FAIL 17 | WARN 0 | SKIP 0 | PASS 114 ]
> 

[pachadotdev]

Hi @apoorvalal !

I could not fix some binomial and Poisson issues (I don't know much Python), so I ended up taking an "average" approach between alpaca, capybara < 2.0.0 and fixest.

I still have to fix some errors regarding APEs. At least this works now for Poisson models and does everything on Armadillo's side.

Btw, have you send pyfixest to @lrberge?