You could create such a matrix with expand.grid(0:1,0:1,0:1,0:1) or, to make it easier with larger numbers of variables, expand.grid(replicate(4, 0:1, simplify=FALSE)), but what about creating the regression formulas instead:
library(tidyverse)
formulas = map(1:4, ~ combn(paste0("x",1:4), .x) %>%
apply(., 2, function(v) paste0("y ~ ", paste(v, collapse=" + ")))) %>%
unlist
formulas
[1] "y ~ x1" "y ~ x2" "y ~ x3"
[4] "y ~ x4" "y ~ x1 + x2" "y ~ x1 + x3"
[7] "y ~ x1 + x4" "y ~ x2 + x3" "y ~ x2 + x4"
[10] "y ~ x3 + x4" "y ~ x1 + x2 + x3" "y ~ x1 + x2 + x4"
[13] "y ~ x1 + x3 + x4" "y ~ x2 + x3 + x4" "y ~ x1 + x2 + x3 + x4"
Then to run the regressions:
models = map(formulas, ~lm(.x, data=df))
There are also some packages with features for doing subsets regression. For example, here.