There is a randomized algorithm that solves this problem in O(n (log n) (log m)) time. It is a Las Vegas algorithm, which means it always give correct results but might take longer than expected. In this case, the probability that it takes much longer than expected is extremely small.

When m = 1, this problem reduces to the problem of finding the median in a read-only array using constant space. That problem does not have a known optimal solution: see "Finding median in read-only memory on integer input, Chan et al."

One odd thing about this reduction of the problem when m = 1 is that this **sub**case is also a **super**case, in that an algorithm for m = 1 can be applied to the m > 1 case. The idea is to just forget that the array rows are sorted and treat the entire storage area as an unsorted array of size n * m. So, for instance, the trivial algorithm for the m = 1 case, in which each element is checked to see if it is the median, takes O(n^{2}) time. Applying it when m > 1 takes O(n^{2}m^{2}) time.

Going back to the m = 1 case, in the comparison model (in which the items of the array can be integers, strings, real numbers, or anything else that can be compared with the inequality operators "<", ">"), the best known deterministic solution that uses space *s* (where *s* is a constant, i.e. in O(1)) has time ϴ(2^{s}s!n^{1 + 1/s}), and it is more complex than the usual algorithms discussed on stackoverflow (though not on https://cstheory.stackexchange.com or https://cs.stackexchange.com). It uses a chained sequence of algorithms A_{s}, A_{s-1}, ..., A_{1}, where A_{s+1} calls A_{s}. You can read it in "Selection from read-only memory and sorting with minimum data movement", by Munro and Raman.

There is a simple randomized algorithm with a smaller run time with high probability. For any constant c, this algorithm runs in time O(n log n) with probability 1 - O(n^{-c}). When the array is the matrix of size n*m that works out to O(n m log (n m)).

This algorithm is very much like quickselect without the rearranging of elements during partitioning.

```
import random
def index_range(needle, haystack):
"""The index range' of a value over an array is a pair
consisting of the number of elements in the array less
than that value and the number of elements in the array
less than or equal to the value.
"""
less = same = 0
for x in haystack:
if x < needle: less += 1
elif x == needle: same += 1
return less, less + same
def median(xs):
"""Finds the median of xs using O(1) extra space. Does not
alter xs.
"""
if not xs: return None
# First, find the minimum and maximum of the array and
# their index ranges:
lo, hi = min(xs), max(xs)
lo_begin, lo_end = index_range(lo, xs)
hi_begin, hi_end = index_range(hi, xs)
# Gradually we will move the lo and hi index ranges closer
# to the median.
mid_idx = len(xs)//2
while True:
print "range size", hi_begin - lo_end
if lo_begin <= mid_idx < lo_end:
return lo
if hi_begin <= mid_idx < hi_end:
return hi
assert hi_begin - lo_end > 0
# Loop over the array, inspecting each item between lo
# and hi. This loops sole purpose is to reservoir sample
# from that set. This makes res a randomly selected
# element from among those strictly between lo and hi in
# xs:
res_size = 0
res = None
for x in xs:
if lo < x < hi:
res_size += 1
if 1 == random.randint(1, res_size):
res = x
assert res is not None
assert hi_begin - lo_end == res_size
# Now find which size of the median res is on and
# continue the search on the smaller region:
res_begin, res_end = index_range(res, xs)
if res_end > mid_idx:
hi, hi_begin, hi_end = res, res_begin, res_end
else:
lo, lo_begin, lo_end = res, res_begin, res_end
```

It works by maintaining upper and lower bounds on the value of the median. It then loops over the array and randomly selects a value between the bounds. That value replaces one of the bounds and the process starts again.

The bounds are accompanied by their index range, a measure of which indexes the bound would appear at if the array were sorted. Once one of the bounds would appear at the index ⌊n/2⌋, it is the median and the algorithm terminates.

When an element is randomly selected in the gap between the bounds, this reduces the gap by 50% in expectation. The algorithm terminates (at the latest) when the gap is 0. We can model this as a series of random independent uniformly distributed variables X_{i} from (0,1) such that Y_{k} = X_{1} * X_{2} * ... * X_{k} where X_{i} is the ratio of the gap that remains after round i. For instance, if after the 10th round the gap between the index ranges of `lo`

and `hi`

is 120, and after the 11th round the gap is 90, then X_{11} = 0.75. The algorithm terminates when Y_{k} < 1/n, because the gap is then less than 1.

Pick a constant positive integer k. Let's bound the probability that Y_{k log2n} >= 1/n using Chernoff bounds. We have Y_{k log2n} = X_{1} * X_{2} * ... X_{k log2n}, so ln Y_{k log2n} = ln X_{1} + ln X_{2} + ... + ln X_{k log2n}. The Chernoff bound then gives Pr(ln X_{1} + ln X_{2} + ... + ln X_{k log2n} >= ln (1/n)) <= min_{t > 0} e^{-t ln (1/n)} (E[e^{t ln X1}] * E[e^{t ln X2}] * ... * E[e^{t ln Xk log2 n}]). After some simplification, the right-hand side is min_{t > 0} n^{t} (E[X_{1}^{t}] * E[X_{2}^{t}] * ... * E[X_{k log2 n}^{t}]). Since this is a minimum and we are looking for an upper bound, we can weaken this by specializing to t = 1. It then simplifies to n^{1-k}, since E[X_{i}] = 1/2.

If we pick, for instance, k = 6, then this bounds the probability that there are 6 log_{2}n rounds or more by n^{-5}. So with probability 1 - O(n^{-5}) the algorithm performs 6 log_{2}n - 1 or fewer rounds. This is what I mean by "with high probability" above.

Since each round inspects every member of the array a constant number of times, each round takes linear time, for a total running time of O(n log n) with high probability. When the array is not just an array but a matrix of size n * m that works out to O(n m log (n m)).

We can do substantially better, however, by taking advantage of the sortedness of the rows. When we were working in a single unsorted array, finding the elements in the gap I referenced above required inspecting each element of the array. In a matrix with sorted rows, the elements in the gap are located in a contiguous segment of each row. Each segment can be identified in O(log m) time using binary search, so they can all be located in O(n log m) time. The reservoir sampling now takes O(n log m) time per iteration of the loop.

The other main work done in the loop is to identify the index range of the element from the gap that was randomly selected. Again, because each row is sorted, the index range for the randomly-chosen element in a row can be determined in O(log m) time. The sums of the index ranges for each row constitute the index range over the whole array, so this part of each loop iteration also takes only O(n log m) time.

By the same argument as above with the Chernoff bound, there are O(log n) iterations with probability at least 1-O(n^{-k}) for any constant k. Thus the whole algorithm takes O(n (log n) (log m)) time with high probability.

```
import bisect
import random
def matrix_index_range(needle, haystack):
"""matrix_index_range calculates the index range of needle
in a haystack that is a matrix (stored in row-major order)
in which each row is sorted"""
n, m = len(haystack), len(haystack[0])
begin = end = 0;
for x in haystack:
begin += bisect.bisect_left(x, needle)
end += bisect.bisect_right(x, needle)
return begin, end
def matrix_median(xs):
print "Starting"
if not xs or not xs[0]: return None
n, m = len(xs), len(xs[0])
lo, hi = xs[0][0], xs[0][m-1]
for x in xs:
lo, hi = min(lo, x[0]), max(hi, x[m-1])
lo_begin, lo_end = matrix_index_range(lo, xs)
hi_begin, hi_end = matrix_index_range(hi, xs)
mid_idx = (n * m) // 2
while True:
print "range size", hi_begin - lo_end
if lo_begin <= mid_idx < lo_end:
return lo
if hi_begin <= mid_idx < hi_end:
return hi
assert hi_begin - lo_end > 0
mid = None
midth = random.randint(0, hi_begin - lo_end - 1)
for x in xs:
gap_begin = bisect.bisect_right(x, lo)
gap_end = bisect.bisect_left(x, hi)
gap_size = gap_end - gap_begin
if midth < gap_size:
mid = x[gap_begin + midth]
break
midth -= gap_size
assert mid is not None
mid_begin, mid_end = matrix_index_range(mid, xs)
assert lo_end <= mid_begin and mid_end <= hi_begin
if mid_end > mid_idx:
hi, hi_begin, hi_end = mid, mid_begin, mid_end
else:
lo, lo_begin, lo_end = mid, mid_begin, mid_end
```

This solution is substantially faster than the first one when m is non-constant.