R - min, max and mean of off-diagonal elements in a matrix
Here the row() and col() helper functions are useful. Using @James A, we can get the upper off-diagonal using this little trick:
> A[row(A) == (col(A) - 1)]
[1] 5 10 15
and the lower off diagonal via this:
> A[row(A) == (col(A) + 1)]
[1] 2 7 12
These can be generalised to give whatever diagonals you want:
> A[row(A) == (col(A) - 2)]
[1] 9 14
and don't require any subsetting.
Then it is a simple matter of calling whatever function you want on these values. E.g.:
> mean(A[row(A) == (col(A) - 1)])
[1] 10
If as per my comment you mean everything but the diagonal, then use
> diag(A) <- NA
> mean(A, na.rm = TRUE)
[1] 8.5
> max(A, na.rm = TRUE)
[1] 15
> # etc. using sum(A, na.rm = TRUE), min(A, na.rm = TRUE), etc..
So this doesn't get lost, Ben Bolker suggests (in the comments) that the above code block can be done more neatly using the row() and col() functions I mentioned above:
mean(A[row(A)!=col(A)])
min(A[row(A)!=col(A)])
max(A[row(A)!=col(A)])
sum(A[row(A)!=col(A)])
which is a nicer solution all round.
Extract opposite-diagonal (not off-diagonal) elements of a matrix
An approach could be
r[(n<-nrow(r))^2-(1:n)*(n-1)]
# [1] 7 2 8
## microbenchmark (matrix(1:1e6,1000))
# Unit: microseconds
# expr min lq mean median uq max neval
# r[(n<-nr... 26.897 39.0075 65.36835 47.309 85.9345 316.97 100
# diag(r[,... 18070.388 18905.3475 20237.09599 19956.615 20423.4695 27798.88 100
# rev(r[ro... 14220.609 21206.7220 21238.59515 22036.275 22599.4490 33252.58 100
Get all diagonal vectors from matrix
A <- matrix(1:16, 4)
# create an indicator for all diagonals in the matrix
d <- row(A) - col(A)
# use split to group on these values
split(A, d)
#
# $`-3`
# [1] 13
#
# $`-2`
# [1] 9 14
#
# $`-1`
# [1] 5 10 15
#
# $`0`
# [1] 1 6 11 16
#
# $`1`
# [1] 2 7 12
#
# $`2`
# [1] 3 8
#
# $`3`
# [1] 4
Find one non-diagonal minimum element from a matrix
For a fully vectorized alternative try
B = (A + diag(Inf(size(diag(A))))); % put Inf on diagonal
[~,ndx] = min(B(:)); % get the linear index of the minimum value
[r,c] = ind2sub(size(A),ndx); % get row, column of corresponding to linear index
Zeroing some entries of matrix parallel to matrix diagonal in R
set.seed(42)
d <- matrix(rnorm(25), nrow=5)
zero.them<-function(x, n){
x2<-x
diag(x)<-0
for(i in 1:n){
diag(x[,(-1:-i)])<-0
diag(x[(-1:-i),])<-0
}
x2[which(x==x2)] <- 0
return(x2)
}
zero.them(d, 1) # for 1 above and below the diagonal
[,1] [,2] [,3] [,4] [,5]
[1,] 1.3709584 -0.10612452 0.0000000 0.000000 0.000000
[2,] -0.5646982 1.51152200 2.2866454 0.000000 0.000000
[3,] 0.0000000 -0.09465904 -1.3888607 -2.656455 0.000000
[4,] 0.0000000 0.00000000 -0.2787888 -2.440467 1.214675
[5,] 0.0000000 0.00000000 0.0000000 1.320113 1.895193
zero.them(d, 2) # for 2 above and below the diagonal
[,1] [,2] [,3] [,4] [,5]
[1,] 1.3709584 -0.10612452 1.3048697 0.0000000 0.0000000
[2,] -0.5646982 1.51152200 2.2866454 -0.2842529 0.0000000
[3,] 0.3631284 -0.09465904 -1.3888607 -2.6564554 -0.1719174
[4,] 0.0000000 2.01842371 -0.2787888 -2.4404669 1.2146747
[5,] 0.0000000 0.00000000 -0.1333213 1.3201133 1.8951935
Find maximum and minimum nondiagonal elements of matrix
The matrix must be at least a 3x3 matrix to have off-diagonal elements. If you have an NxN matrix, then (using r for row and c for column) the diagonals have r == c and r == N - c - 1 (aka r + c == N - 1). You can simply ignore elements for which either of those conditions is true.
You must place the other elements into one of the four quadrants aligned with the diagonals. The triangle at the top is the N (north) quadrant; at the bottom is the S (south) quadrant; at the left is the W (west) quadrant; at the right is the E (east) quadrant.
N=5 N=4
r\c 0 1 2 3 4 r\c 0 1 2 3
0 . N N N . 0 . N N .
1 W . N . E 1 W . . E
2 W W . E E 2 W . . E
3 W . S . E 3 . S S .
4 . S S S .
The dots represent the diagonals. In the 5x5 matrix, the W and S elements satisfy r > c; the N and E elements satisfy r < c (and the elements on the diagonal satisfy r == c). Similarly, the W and N elements satisfy r + c < N - 1 while the S and E elements satisfy r + c > N -1 (and the elements on the diagonal satisfy r + c == N - 1). The 4x4 matrix also satisfies these conditions.
You can then use these conditions to classify a given (row, column) pair into one of the four quadrants using an encoding scheme similar to that suggested by Dèja vu in a comment, but adjusted to use different conditions:
int quadrant = ((r < c) ? 0 : 1) + ((r + c < N - 1) ? 0 : 2);
When quadrant is 0, it corresponds to the N quadrant; when it is 1, to the W quadrant; when it is 2, to the E quadrant; when it is 3, to the S quadrant. Note that this condition assumes that the diagonal elements are not tested. You might write assert(r != c && r + c != N - 1); to ensure that it is not invoked for an element on the diagonals.
You can then compare the element mat[r][c] with min[quadrant] and max[quadrant], which is the scheme outlined by kaylum in a comment.
Possible code (matminmax37.c)
/* SO 7085-8579 */
#include <assert.h>
#include <stdio.h>
static void print_matrix(const char *tag, size_t n, size_t m, const int matrix[n][m]);
/*
** When quadrant is 0, it corresponds to the N quadrant; when it is 1,
** to the W quadrant; 2, to the E quadrant; when it is 3, to the S
** quadrant. Note that this condition assumes that the diagonal
** elements are not tested; the assertion ensures that.
*/
static inline size_t quadrant(size_t r, size_t c, size_t N)
{
assert(r != c && r + c != N - 1);
return ((r < c) ? 0 : 1) + ((r + c < N - 1) ? 0 : 2);
}
static void matminmax(size_t n, const int matrix[n][n], int min[4], int max[4])
{
min[0] = max[0] = matrix[0][1];
min[1] = max[1] = matrix[1][0];
min[2] = max[2] = matrix[1][n-1];
min[3] = max[3] = matrix[n-1][1];
for (size_t r = 0; r < n; r++)
{
for (size_t c = 0; c < n; c++)
{
if (r == c || r + c == n - 1)
continue;
size_t q = quadrant(r, c, n);
if (matrix[r][c] > max[q])
max[q] = matrix[r][c];
if (matrix[r][c] < min[q])
min[q] = matrix[r][c];
}
}
}
static void print_minmax(int min[4], int max[4])
{
static const char quadname[] = "NWES";
putchar('\n');
for (size_t i = 0; i < 4; i++)
printf("%c quadrant: min = %4d, max = %4d\n", quadname[i], min[i], max[i]);
putchar('\n');
}
int main(void)
{
/* Created by: gen_matrix -C -H 999 -L -999 -c 12 -r 12 -i -n mtx1 -w 4 */
/* Random seed: 0xDB2DC719 */
const int mtx1[12][12] =
{
{ 688, 351, -309, -491, -200, 958, -699, -284, -638, 696, -427, 43, },
{ -621, 242, 913, -247, -293, 878, 538, 368, -861, 488, -568, 394, },
{ -646, -813, -943, -417, -697, 551, 325, 11, -234, 18, 348, -229, },
{ 436, 529, -994, 457, -270, 369, 318, -455, -918, 444, -821, 266, },
{ 544, 24, 398, -589, -779, -2, 932, 810, -915, -591, 877, 865, },
{ 19, 879, -528, -483, 961, 478, -426, 528, 522, -403, -859, -63, },
{ 830, 103, -662, -894, -29, 875, -610, 244, 592, -28, -206, 538, },
{ 231, -707, -303, -462, 683, 200, -3, -707, 480, -2, -920, -743, },
{ 981, -427, -793, -577, -296, 127, 823, 937, -98, -405, 871, -221, },
{ -139, 841, 299, -415, -773, 403, 92, -474, -384, -219, -636, -12, },
{ -89, -13, -679, -998, -567, -826, 705, -530, 658, 134, 871, 875, },
{ -475, -148, -275, 754, 399, -946, -78, -469, 655, -45, 81, 921, },
};
enum { MTX1_SIZE = sizeof(mtx1) / sizeof(mtx1[0]) };
int min[4];
int max[4];
matminmax(MTX1_SIZE, mtx1, min, max);
print_matrix("Matrix 1", MTX1_SIZE, MTX1_SIZE, mtx1);
print_minmax(min, max);
/* Created by: gen_matrix -C -H 999 -L -999 -c 13 -r 13 -i -n mtx2 -w 4 */
/* Random seed: 0xD2A82AC2 */
const int mtx2[13][13] =
{
{ 208, -676, 374, 587, -89, -485, -754, 286, -295, -826, -511, -797, 858, },
{ 991, 423, 404, 477, -449, 442, -860, 629, 437, -606, 974, 522, 885, },
{ -54, -815, -124, 27, 68, -224, -95, 430, -244, -941, 857, -843, -306, },
{ 92, 941, -613, 435, 222, -966, -1, 292, -577, 597, 238, 527, -984, },
{ 978, 661, 315, -347, -747, 242, 711, 15, -922, -623, -533, 794, 65, },
{ 459, -195, 928, -325, -270, -703, 64, -18, -219, 92, 831, -657, 945, },
{ 250, -102, -861, -30, -603, -921, -229, 914, 164, 273, -133, 915, 565, },
{ 855, -607, 800, -137, 635, -216, -990, -432, -986, -650, 850, 456, 671, },
{ -393, -159, -685, 537, 598, -680, 241, -799, -821, -750, -559, -164, 103, },
{ 318, 679, 499, 980, 186, 841, -139, 878, -187, 818, -430, 904, 797, },
{ -496, 31, -606, 25, -244, 26, -558, -307, -656, 885, 327, -961, -742, },
{ 558, 592, -562, 148, -543, -995, 534, -519, 257, 784, 982, 508, 470, },
{ -37, 289, -94, -513, -115, -51, 638, 633, 64, -701, -955, -575, 987, },
};
enum { MTX2_SIZE = sizeof(mtx2) / sizeof(mtx2[0]) };
matminmax(MTX2_SIZE, mtx2, min, max);
print_matrix("Matrix 2", MTX2_SIZE, MTX2_SIZE, mtx2);
print_minmax(min, max);
/* Created by: gen_matrix -C -H 999 -L -999 -c 3 -r 3 -i -n mtx3 -w 4 */
/* Random seed: 0x238E2E67 */
const int mtx3[3][3] =
{
{ -703, 420, 896, },
{ 697, -581, -38, },
{ 829, 878, 722, },
};
enum { MTX3_SIZE = sizeof(mtx3) / sizeof(mtx3[0]) };
matminmax(MTX3_SIZE, mtx3, min, max);
print_matrix("Matrix 3", MTX3_SIZE, MTX3_SIZE, mtx3);
print_minmax(min, max);
/* Created by: gen_matrix -C -H 999 -L -999 -c 4 -r 4 -i -n mtx4 -w 4 */
/* Random seed: 0x9F26F797 */
const int mtx4[4][4] =
{
{ -726, 135, 216, 291, },
{ -628, -901, -594, 313, },
{ -315, -696, 274, 337, },
{ 499, 585, -148, -390, },
};
enum { MTX4_SIZE = sizeof(mtx4) / sizeof(mtx4[0]) };
matminmax(MTX4_SIZE, mtx4, min, max);
print_matrix("Matrix 4", MTX4_SIZE, MTX4_SIZE, mtx4);
print_minmax(min, max);
/* Created by: gen_matrix -C -H 999 -L -999 -c 5 -r 5 -i -n mtx5 -w 4 */
/* Random seed: 0x10185C02 */
const int mtx5[5][5] =
{
{ -73, -575, -606, 445, -714, },
{ 995, -839, -773, 240, 819, },
{ 132, 581, 956, -495, -914, },
{ -10, -748, 328, -611, -475, },
{ 249, 249, -839, 152, -408, },
};
enum { MTX5_SIZE = sizeof(mtx5) / sizeof(mtx5[0]) };
matminmax(MTX5_SIZE, mtx5, min, max);
print_matrix("Matrix 5", MTX5_SIZE, MTX5_SIZE, mtx5);
print_minmax(min, max);
/* Created by: gen_matrix -C -H 999 -L -999 -c 6 -r 6 -i -n mtx6 -w 4 */
/* Random seed: 0x875034FD */
const int mtx6[6][6] =
{
{ -814, -712, -938, -333, 123, -626, },
{ 621, -154, 835, 640, 575, -287, },
{ -720, -288, -563, 291, 182, -542, },
{ 132, -954, -404, -859, -796, 212, },
{ -506, 237, -828, -37, 431, 399, },
{ -943, 151, 567, -414, -902, 959, },
};
enum { MTX6_SIZE = sizeof(mtx6) / sizeof(mtx6[0]) };
matminmax(MTX6_SIZE, mtx6, min, max);
print_matrix("Matrix 6", MTX6_SIZE, MTX6_SIZE, mtx6);
print_minmax(min, max);
/* Created by: gen_matrix -C -H 999 -L -999 -c 7 -r 7 -i -n mtx7 -w 4 */
/* Random seed: 0x2574C5FF */
const int mtx7[7][7] =
{
{ 386, -560, 959, -485, 360, 278, -474, },
{ 440, -360, 617, 214, -11, 300, 590, },
{ -285, 83, -568, 989, -117, -766, 411, },
{ 944, 598, -568, 414, -596, -882, 640, },
{ -544, -849, 167, 890, 474, -361, 52, },
{ 531, 544, -488, 569, 585, -816, 486, },
{ -542, -965, -758, 590, -314, 401, 964, },
};
enum { MTX7_SIZE = sizeof(mtx7) / sizeof(mtx7[0]) };
matminmax(MTX7_SIZE, mtx7, min, max);
print_matrix("Matrix 7", MTX7_SIZE, MTX7_SIZE, mtx7);
print_minmax(min, max);
return 0;
}
static int max_field_width(size_t n, const int matrix[n][n])
{
int min_val = matrix[0][0];
int max_val = matrix[0][0];
for (size_t i = 0; i < n; i++)
{
for (size_t j = 0; j < n; j++)
{
if (matrix[i][j] < min_val)
min_val = matrix[i][j];
if (matrix[i][j] > max_val)
max_val = matrix[i][j];
}
}
int fld_width = snprintf(0, 0, "%d", max_val);
if (min_val < 0)
{
int min_width = snprintf(0, 0, "%d", min_val);
if (min_width > fld_width)
fld_width = min_width;
}
return fld_width;
}
static void print_matrix(const char *tag, size_t n, size_t m, const int matrix[n][m])
{
printf("%s (%zux%zu):\n", tag, n, m);
int w = max_field_width(n, matrix) + 1;
for (size_t i = 0; i < n; i++)
{
for (size_t j = 0; j < m; j++)
{
printf("%*d", w, matrix[i][j]);
}
putchar('\n');
}
}
Output
Matrix 1 (12x12):
688 351 -309 -491 -200 958 -699 -284 -638 696 -427 43
-621 242 913 -247 -293 878 538 368 -861 488 -568 394
-646 -813 -943 -417 -697 551 325 11 -234 18 348 -229
436 529 -994 457 -270 369 318 -455 -918 444 -821 266
544 24 398 -589 -779 -2 932 810 -915 -591 877 865
19 879 -528 -483 961 478 -426 528 522 -403 -859 -63
830 103 -662 -894 -29 875 -610 244 592 -28 -206 538
231 -707 -303 -462 683 200 -3 -707 480 -2 -920 -743
981 -427 -793 -577 -296 127 823 937 -98 -405 871 -221
-139 841 299 -415 -773 403 92 -474 -384 -219 -636 -12
-89 -13 -679 -998 -567 -826 705 -530 658 134 871 875
-475 -148 -275 754 399 -946 -78 -469 655 -45 81 921
N quadrant: min = -861, max = 958
W quadrant: min = -994, max = 981
E quadrant: min = -920, max = 877
S quadrant: min = -998, max = 937
Matrix 2 (13x13):
208 -676 374 587 -89 -485 -754 286 -295 -826 -511 -797 858
991 423 404 477 -449 442 -860 629 437 -606 974 522 885
-54 -815 -124 27 68 -224 -95 430 -244 -941 857 -843 -306
92 941 -613 435 222 -966 -1 292 -577 597 238 527 -984
978 661 315 -347 -747 242 711 15 -922 -623 -533 794 65
459 -195 928 -325 -270 -703 64 -18 -219 92 831 -657 945
250 -102 -861 -30 -603 -921 -229 914 164 273 -133 915 565
855 -607 800 -137 635 -216 -990 -432 -986 -650 850 456 671
-393 -159 -685 537 598 -680 241 -799 -821 -750 -559 -164 103
318 679 499 980 186 841 -139 878 -187 818 -430 904 797
-496 31 -606 25 -244 26 -558 -307 -656 885 327 -961 -742
558 592 -562 148 -543 -995 534 -519 257 784 982 508 470
-37 289 -94 -513 -115 -51 638 633 64 -701 -955 -575 987
N quadrant: min = -966, max = 974
W quadrant: min = -921, max = 991
E quadrant: min = -986, max = 945
S quadrant: min = -995, max = 982
Matrix 3 (3x3):
-703 420 896
697 -581 -38
829 878 722
N quadrant: min = 420, max = 420
W quadrant: min = 697, max = 697
E quadrant: min = -38, max = -38
S quadrant: min = 878, max = 878
Matrix 4 (4x4):
-726 135 216 291
-628 -901 -594 313
-315 -696 274 337
499 585 -148 -390
N quadrant: min = 135, max = 216
W quadrant: min = -628, max = -315
E quadrant: min = 313, max = 337
S quadrant: min = -148, max = 585
Matrix 5 (5x5):
-73 -575 -606 445 -714
995 -839 -773 240 819
132 581 956 -495 -914
-10 -748 328 -611 -475
249 249 -839 152 -408
N quadrant: min = -773, max = 445
W quadrant: min = -10, max = 995
E quadrant: min = -914, max = 819
S quadrant: min = -839, max = 328
Matrix 6 (6x6):
-814 -712 -938 -333 123 -626
621 -154 835 640 575 -287
-720 -288 -563 291 182 -542
132 -954 -404 -859 -796 212
-506 237 -828 -37 431 399
-943 151 567 -414 -902 959
N quadrant: min = -938, max = 835
W quadrant: min = -954, max = 621
E quadrant: min = -796, max = 399
S quadrant: min = -902, max = 567
Matrix 7 (7x7):
386 -560 959 -485 360 278 -474
440 -360 617 214 -11 300 590
-285 83 -568 989 -117 -766 411
944 598 -568 414 -596 -882 640
-544 -849 167 890 474 -361 52
531 544 -488 569 585 -816 486
-542 -965 -758 590 -314 401 964
N quadrant: min = -560, max = 989
W quadrant: min = -849, max = 944
E quadrant: min = -882, max = 640
S quadrant: min = -965, max = 890
Extract diagonals from a distance matrix in R
One work around is to convert the dist object to matrix and then extract elements where row index is one larger than the column index:
mat = as.matrix(dist(mymatrix))
mat[row(mat) == col(mat) + 1]
# [1] 2.828427 3.000000 2.828427
How to find offset diagonal of a matrix?
1) The i-th offset reverse diagonal of square matrix m:
off.rev.diag <- function(m, i = 0) m[ (row(m) + col(m) - 1) %% ncol(m) == i ]
For example:
> m <- matrix(1:25, 5); m
[,1] [,2] [,3] [,4] [,5]
[1,] 1 6 11 16 21
[2,] 2 7 12 17 22
[3,] 3 8 13 18 23
[4,] 4 9 14 19 24
[5,] 5 10 15 20 25
> off.rev.diag(m, 1)
[1] 1 10 14 18 22
> off.rev.diag(m, 2)
[1] 2 6 15 19 23
2) We can also write a replacement function:
"off.rev.diag<-" <- function(m, i = 0, value) {
m.ix <- matrix(seq_along(m), nrow(m))
replace(m, off.rev.diag(m.ix, i), value)
}
For example,
> off.rev.diag(m, 1) <- -(1:5)
> m
[,1] [,2] [,3] [,4] [,5]
[1,] -1 6 11 16 21
[2,] 2 7 12 17 -5
[3,] 3 8 13 -4 23
[4,] 4 9 -3 19 24
[5,] 5 -2 15 20 25
Algorithm to maximize the smallest diagonal element of a matrix
Suppose you wanted to know whether there was a better solution to your problem than 3.
Change your matrix to have a 1 for every element that is strictly greater than 3:
4 3 2 1 1 0 0 0
1 4 3 2 0 1 0 0
2.5 3.5 4.5 1.5 -> 0 1 1 0
2 1 4 3 0 0 1 0
Your problem can be interpreted as trying to find a perfect matching in the bipartite graph which has this binary matrix as its biadjacency graph.
In this case, it is easy to see that there is no way of improving your result because there is no way of reordering rows to make the diagonal entry in the last column greater than 3.
For a larger matrix, there are efficient algorithms to determine maximal matchings in bipartite graphs.
This suggests an algorithm:
- Use bisection to find the largest value for which the generated graph has a perfect matching
- The assignment corresponding to the perfect matching with the largest value will be equal to the best permutation of rows
EDIT
This Python code illustrates how to use the networkx library to determine whether the graph has a perfect matching for a particular cutoff value.
import networkx as nx
A = [[4,3,2,1],
[1,4,3,2],
[2,1,4,3],
[2.5,3.5,4.5,1.5]]
cutoff = 3
G=nx.DiGraph()
for i,row in enumerate(A):
G.add_edge('start','row'+str(i),capacity=1.0)
G.add_edge('col'+str(i),'end',capacity=1.0)
for j,e in enumerate(row):
if e>cutoff:
G.add_edge('row'+str(i),'col'+str(j),capacity=1.0)
if nx.max_flow(G,'start','end')<len(A):
print 'No perfect matching'
else:
print 'Has a perfect matching'
For a random matrix of size 1000*1000 it takes about 1 second on my computer.
Related Topics
Creating Igraph with Isolated Nodes
Plot Event Sequences/Event Sequences Clustering
How to Make Shinyapp to Use Environmental Variables When Deployed on the Web
Specify Function Parameters in Do.Call
How to Scale the Size of Line and Point Separately in Ggplot2
How to Make Join Operations in Dplyr Silent
Independently Move 2 Legends Ggplot2 on a Map
Pass String as Name of Attached Data Column Name
Control Number Formatting in Shiny's Implementation of Datatable
Gradient Breaks in a Ggplot Stat_Bin2D Plot
How to Read the Files in a Directory in Sorted Order Using R
How to Bookmark and Restore Dynamically Added Modules
Importing S3 Method from Another Package
How to Load CSV File into Sparkr on Rstudio
Adding a Ranking Column to a Dataframe
Two Y Axis in Highcharter in R
Error in Bind_Rows_(X, .Id):Column Can't Be Converted from Factor to Numeric