I'm not completely satisfied with this answer, but I think it's better than no answer. I'll attempt to answer the question posed in the comments: in a 100-by-100 matrix of independent Bernoulli trials ($M$), what is the probability of finding a 5-by-5 submatrix, $X$, composed entirely of successes? I use a 0.5 probability of the Bernoulli trials returning 1, but the methodology should work for any probability. The code below is in R.
First, I used a couple Rcpp functions. The first is a fast Rcpp function saved in "C:/temp/findSquare.cpp". It counts how many matrices contain a square submatrix of size n composed entirely of 1s. Each column of the input matrix x holds the values of an nr-by-nc matrix.
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
int findSquare(IntegerMatrix& x, const int& nr, const int& nc, const int& n) {
const int iter = x.nrow();
int out = 0;
for (int k = 0; k < iter; k++) {
int jj = nc;
int cnt = 0;
int found = 0;
for (int j = 0; j < nc; j++) {
cnt = (cnt + 1)*x(k, j);
x(k, j) = (cnt >= n) + 1;
}
for (int i = 1; i < nr; i++) {
int cnt = 0;
for (int j = 0; j < nc; j++) {
cnt = cnt*x(k, jj) + x(k, jj);
x(k, j) = x(k, j)*(cnt >= n) + 1;
found += x(k, j) > n;
jj++;
}
}
out += found > 0;
}
return(out);
}
The second is for matrix multiplication since matrix multiplication in base R is so slow on my setup. It is saved in "C:/temp/eigenMult.cpp"
// [[Rcpp::depends(RcppArmadillo, RcppEigen)]]
#include <RcppArmadillo.h>
#include <RcppEigen.h>
// [[Rcpp::export]]
SEXP eigenMapMatSquare(const Eigen::Map<Eigen::MatrixXd> X){
Eigen::MatrixXd X2 = X * X;
return Rcpp::wrap(X2);
}
// [[Rcpp::export]]
SEXP eigenMapMatMult(const Eigen::Map<Eigen::MatrixXd> A, Eigen::Map<Eigen::MatrixXd> B){
Eigen::MatrixXd C = A * B;
return Rcpp::wrap(C);
}
Finally, an R function for raising a matrix to a power using "eigenMult.cpp":
matpow <- function(x, pow) {
if (pow == 0L) {
1L
} else if(pow %% 2L) {
nextpow <- (pow - 1L)/2L
if (nextpow) {
eigenMapMatMult(x, Recall(eigenMapMatSquare(x), nextpow))
} else x
} else {
Recall(eigenMapMatSquare(x), pow/2L)
}
}
We can get a credible interval via Monte Carlo simulation:
Rcpp::sourceCpp("C:/temp/findSquare.cpp")
# Monte Carlo - take 100 sets of 10k samples
system.time(counts <- sum(replicate(100, findSquare(matrix(sample(0:1, 1e8, TRUE), 1e4), 100L, 100L, 5L))))
#> user system elapsed
#> 641.66 27.36 669.52
# 90% credible interval
(pSimBounds <- qbeta(c(0.05, 0.95), sum(counts) + 0.5, 1e6 - sum(counts) + 0.5))
#> [1] 0.0002212588 0.0002728769
These bounds are fairly tight, but they took more than 11 minutes to compute. A simple Markov chain will give quick and easy bounds. Instead of a 100-by-100 matrix, look for $X$ in a 100-by-5 matrix. The probability of getting 1 for all five values in a row is $2^{-5}$. The transition matrix is then:
m <- cbind(matrix(c(rep(1 - 1/32, 5), 0), 6), diag(1/32, 6)[,-6])
m[6, 6] <- 1
Get the probability of finding $X$:
Rcpp::sourceCpp("C:/temp/eigenMult.cpp")
(p1 <- matpow(m, 100)[1, 6])
#> [1] 2.772544e-06
The full matrix has 20 independent sets of 5 columns. A lower bound could be the probability of finding $X$ only in those 20 independent sets. For an upper bound, there are 96 sets of 5 columns in the full matrix, but they aren't independent: if the $X$ is not in columns 1 through 5, it is less likely to be in columns 2 through 6, etc. Treating the 96 sets as independent will over-estimate the probability of finding $X$ in $M$:
(pLower <- 1 - (1 - p1)^20)
#> [1] 5.544942e-05
(pUpper <- 1 - (1 - p1)^96)
#> [1] 0.0002661292
The lower bound is not very impressive, but the upper bound is tighter than what was calculated from 1M samples.
For the final estimate, let's examine what happens if we start with 5 columns and repeatedly add one more column. For six columns, the Markov chain needs to track what is happening in columns 1:5 as well as 2:6. Each set can be in state 0 through 4 (the number of consecutive rows of 5 1s), resulting in 25 states plus 1 absorbing state. As an example, the sequence 0,1,1,1,1,1 in a row will transition columns 1:5 to state 0 and increment the state of columns 2:6 by 1 (unless the system is in the absorbing state).
If we start with 5 columns and add 1 column at a time, we can calculate the probability of finding $X$ in the last 5 columns conditioned on $X$ not existing in the matrix without the last column. As more columns are added, this probability quickly converges, which allows us to estimate the probability of finding $X$ in $M$. The size of the transition matrix grows quickly as more columns are added, so I stop at 11.
# more involved estimate from multiple Markov chains
library(data.table)
library(Rfast) # for eachrow and rowMaxs functions
fTrans <- function(dt, mTrans, transk, n, nn) {
# iteratively build out state changes
nmTo <- paste0("to", 1:nn)
nmFrom <- colnames(dt)
nmAll <- c(nmFrom, nmTo)
dtStates <- setnames(as.data.table(rep(list(integer(0)), nn*2L + 1L)), c(nmAll, "k"))
while (nrow(dt)) {
ldt <- vector("list", nrow(dt) + 1L)
ldt[1:nrow(dt)] <- lapply(
1:nrow(dt),
function(i) {
m2 <- eachrow(mTrans, as.numeric(dt[i])) + mTrans
m2[rowMaxs(m2, TRUE) == n, ] <- n
cbind(
setnames(dt[rep(i, nrow(mTrans))], nmFrom),
setnames(as.data.table(m2), nmTo)
)[, k := transk][, .(k = sum(k)), nmAll]
}
)
ldt[[nrow(dt) + 1L]] <- dtStates
dtStates <- rbindlist(ldt)
dt <- rbindlist(
list(
unique(dtStates[, 1:nn])[, k := 2L],
unique(dtStates[, (nn + 1L):(2L*nn)])[, k := 1L]
), FALSE
)[
, .(k = sum(k)), nmFrom
][
k == 1L & from1 < n
][
, k := NULL
]
}
dtStates
}
pGridInArray <- function(n, dims, p = 0.5) {
# main function
if (dims[1] < dims[2]) dims <- dims[2:1]
n1 <- n - 1L
nn <- dims[2] - n1 # number of column sets
nn1 <- nn - 1L
nn2 <- 2^nn1
# get all possible row values that will not lead to a transition back to the
# initial state
m <- as.matrix(expand.grid(rep(list(0:1), nn1)))
m <- unique(
do.call(
rbind,
lapply(
1:nn,
function(i) {
mm <- matrix(1L, nn2, dims[2])
if (i > 1L) mm[, 1:(i - 1L)] <- m[, 1:(i - 1L)]
if (i + n <= dims[2]) mm[, (i + n):dims[2]] <- m[, i:nn1]
mm
}
)
)
)
# for each row in m, which column sets would it increment?
# aggregate unique counts to feed the creation of a transition matrix
dtM <- as.data.table(
m %*% (diag(0L, dims[2], nn) + rep(c(rep(1L, n), integer(nn)), nn)[1:(dims[2]*nn)]) == n
)[
, .(k = .N), eval(paste0("V", 1:nn))
]
mTrans <- as.matrix(dtM[, 1:nn])
class(mTrans) <- "numeric"
# build a data.table listing all the state changes that don't end in the
# initial state, and the number of ways the state change can occur
dtTrans <- setorder(
fTrans(
setnames(
as.data.table(as.list(integer(nn))),
paste0("from", 1:nn)
),
mTrans, dtM$k, n, nn
)
)
dtTrans[
, `:=`(
state1 = rleidv(dtTrans[, 1:nn]),
state2 = setorder(
setorder(
dtTrans[, (nn + 1L):(2L*nn)][, idx := .I]
)[
, state2 := rleidv(.SD) + 1L, .SDcols = 1:nn
], idx
)$state2
)
]
# build the transition matrix, calculate dims[1] transitions, and get the
# probability of ending in the absorbing state (when the desired submatrix has
# occurred)
nStates <- max(dtTrans$state2)
m <- matrix(0, nStates, nStates)
m[as.matrix(dtTrans[, state1:state2])] <- dtTrans$k
m[nStates, nStates] <- 2^dims[2] # absorbing state
m[, 1] <- 2^dims[2] - sum(dtM$k)
m[nStates, 1] <- 0
matpow(m*p^dims[2], dims[1])[1, nStates]
}
# test the pGridInArray function
pGridInArray(2L, c(5L, 3L))
#> [1] 0.3262634
findSquare(as.matrix(expand.grid(rep(list(0:1), 15))), 5L, 3L, 2L)/2^15
#> [1] 0.3262634
pGridInArray(3L, c(5L, 4L))
#> [1] 0.009971619
findSquare(as.matrix(expand.grid(rep(list(0:1), 20))), 5L, 4L, 3L)/2^20
#> [1] 0.009971619
pGridInArray(5L, c(100L, 5L))
#> [1] 2.772544e-06
p1
#> [1] 2.772544e-06
# 5-by-5 filled submatrix in 100 rows x 5-11 columns
p5 <- vapply(5:11, function(i) pGridInArray(5L, c(100L, i)), 0)
dt <- data.table(
columns = 5:11,
prob = p5,
# probability of finding the submatrix in the last 5 columns
# given that it is not anywhere else
marg = c(NA, diff(p5)/(1 - first(p5, -1)))
)[
# project the probability of finding the submatrix in 100 columns
# using the last "marginal" probability
, proj := 1 - (1 - prob)*(1 - marg)^(100 - columns)
]
dt
#> columns prob marg proj
#> 1 5 2.772544e-06 NA NA
#> 2 6 5.454110e-06 2.681573e-06 0.0002574892
#> 3 7 8.135459e-06 2.681364e-06 0.0002574695
#> 4 8 1.081678e-05 2.681347e-06 0.0002574679
#> 5 9 1.349810e-05 2.681345e-06 0.0002574678
#> 6 10 1.617941e-05 2.681346e-06 0.0002574678
#> 7 11 1.886071e-05 2.681346e-06 0.0002574678
# show that the estimate is converging
diff(last(dt$proj, -1))
#> [1] -1.969634e-08 -1.576093e-09 -1.256927e-10 2.543288e-11 -5.694340e-13
# final estimate
last(dt$proj)
#> [1] 0.0002574678
Unfortunately, I don't have a method to place error bounds on the final estimate, nor can I prove that the estimate is strictly converging with each iteration.
