179 lines
4.7 KiB
R
179 lines
4.7 KiB
R
# ---------------------------------------------------------------
|
||
# Strassen matrix multiplication – R version (Corrected)
|
||
# ---------------------------------------------------------------
|
||
|
||
# Helper function to pad matrix to next power of 2
|
||
pad_to_power_of_2 <- function(mat) {
|
||
n <- nrow(mat)
|
||
m <- ncol(mat)
|
||
max_dim <- max(n, m)
|
||
new_dim <- 2^ceiling(log2(max_dim))
|
||
|
||
if (new_dim == n && new_dim == m) {
|
||
return(mat)
|
||
}
|
||
|
||
padded <- matrix(0, nrow = new_dim, ncol = new_dim)
|
||
padded[1:n, 1:m] <- mat
|
||
return(padded)
|
||
}
|
||
|
||
# Helper to split a matrix into quadrants
|
||
split_matrix <- function(X) {
|
||
n <- nrow(X)
|
||
mid <- n %/% 2
|
||
|
||
list(
|
||
X11 = X[1:mid, 1:mid, drop = FALSE],
|
||
X12 = X[1:mid, (mid + 1):n, drop = FALSE],
|
||
X21 = X[(mid + 1):n, 1:mid, drop = FALSE],
|
||
X22 = X[(mid + 1):n, (mid + 1):n, drop = FALSE]
|
||
)
|
||
}
|
||
|
||
# Join four submatrices into one matrix
|
||
join_matrices <- function(C11, C12, C21, C22) {
|
||
top <- cbind(C11, C12)
|
||
bottom <- cbind(C21, C22)
|
||
rbind(top, bottom)
|
||
}
|
||
|
||
# Strassen Matrix Multiplication
|
||
strassen <- function(A, B) {
|
||
# Check if matrices are conformable for multiplication
|
||
if (ncol(A) != nrow(B)) {
|
||
stop("Matrices cannot be multiplied: incompatible dimensions")
|
||
}
|
||
|
||
# Pad matrices to make them square and powers of 2
|
||
orig_rows <- nrow(A)
|
||
orig_cols <- ncol(B)
|
||
|
||
max_dim <- max(nrow(A), ncol(A), nrow(B), ncol(B))
|
||
padded_dim <- 2^ceiling(log2(max_dim))
|
||
|
||
if (nrow(A) != padded_dim || ncol(A) != padded_dim) {
|
||
A_padded <- matrix(0, nrow = padded_dim, ncol = padded_dim)
|
||
A_padded[1:nrow(A), 1:ncol(A)] <- A
|
||
A <- A_padded
|
||
}
|
||
|
||
if (nrow(B) != padded_dim || ncol(B) != padded_dim) {
|
||
B_padded <- matrix(0, nrow = padded_dim, ncol = padded_dim)
|
||
B_padded[1:nrow(B), 1:ncol(B)] <- B
|
||
B <- B_padded
|
||
}
|
||
|
||
# Recursive helper function
|
||
strassen_recursive <- function(X, Y) {
|
||
n <- nrow(X)
|
||
|
||
# Base case: 1x1 matrix
|
||
if (n == 1) {
|
||
return(matrix(X[1, 1] * Y[1, 1], nrow = 1, ncol = 1))
|
||
}
|
||
|
||
# Split matrices into quadrants
|
||
X_parts <- split_matrix(X)
|
||
Y_parts <- split_matrix(Y)
|
||
|
||
x11 <- X_parts$X11
|
||
x12 <- X_parts$X12
|
||
x21 <- X_parts$X21
|
||
x22 <- X_parts$X22
|
||
|
||
y11 <- Y_parts$X11
|
||
y12 <- Y_parts$X12
|
||
y21 <- Y_parts$X21
|
||
y22 <- Y_parts$X22
|
||
|
||
# Compute the seven products recursively
|
||
p1 <- strassen_recursive(x12 - x22, y21 + y22)
|
||
p2 <- strassen_recursive(x11 + x22, y11 + y22)
|
||
p3 <- strassen_recursive(x11 - x21, y11 + y12)
|
||
p4 <- strassen_recursive(x11 + x12, y22)
|
||
p5 <- strassen_recursive(x11, y12 - y22)
|
||
p6 <- strassen_recursive(x22, y21 - y11)
|
||
p7 <- strassen_recursive(x21 + x22, y11)
|
||
|
||
# Combine results into quadrants
|
||
c11 <- p1 + p2 - p4 + p6
|
||
c12 <- p4 + p5
|
||
c21 <- p6 + p7
|
||
c22 <- p2 - p3 + p5 - p7
|
||
|
||
# Join and return final matrix
|
||
join_matrices(c11, c12, c21, c22)
|
||
}
|
||
|
||
# Perform Strassen multiplication
|
||
result_padded <- strassen_recursive(A, B)
|
||
|
||
# Extract the relevant portion
|
||
result <- result_padded[1:orig_rows, 1:orig_cols, drop = FALSE]
|
||
|
||
return(result)
|
||
}
|
||
|
||
# ---------------------------------------------------------------
|
||
# Helper that mimics the Julia `intprint`
|
||
# ---------------------------------------------------------------
|
||
int_print <- function(title, mat) {
|
||
# Round to 8 decimal places and convert near-integers to integers
|
||
rounded_mat <- round(mat, digits = 8)
|
||
cat(title, "\n")
|
||
print(rounded_mat)
|
||
cat("\n")
|
||
}
|
||
|
||
# ---------------------------------------------------------------
|
||
# Test data
|
||
# ---------------------------------------------------------------
|
||
varA <- matrix(c(1, 2, 3, 4), nrow = 2, byrow = TRUE)
|
||
varB <- matrix(c(5, 6, 7, 8), nrow = 2, byrow = TRUE)
|
||
|
||
varC <- matrix(c(
|
||
1, 1, 1, 1,
|
||
2, 4, 8, 16,
|
||
3, 9, 27, 81,
|
||
4, 16, 64, 256
|
||
), nrow = 4, byrow = TRUE)
|
||
|
||
varD <- matrix(c(
|
||
4, -3, 4/3, -1/4,
|
||
-13/3, 19/4, -7/3, 11/24,
|
||
3/2, -2, 7/6, -1/4,
|
||
-1/6, 1/4, -1/6, 1/24
|
||
), nrow = 4, byrow = TRUE)
|
||
|
||
varE <- matrix(c(
|
||
1, 2, 3, 4,
|
||
5, 6, 7, 8,
|
||
9, 10, 11, 12,
|
||
13, 14, 15, 16
|
||
), nrow = 4, byrow = TRUE)
|
||
|
||
varF <- diag(4)
|
||
|
||
r <- sqrt(2)/2
|
||
R <- matrix(c(r, r, -r, r), nrow = 2, byrow = TRUE)
|
||
|
||
# ---------------------------------------------------------------
|
||
# Run the examples – compare ordinary (%*%) and Strassen results
|
||
# ---------------------------------------------------------------
|
||
|
||
# Example 1
|
||
int_print("Regular multiply:", t(varA) %*% t(varB))
|
||
int_print("Strassen multiply:", strassen(t(varA), t(varB)))
|
||
|
||
# Example 2
|
||
int_print("Regular multiply:", varC %*% varD)
|
||
int_print("Strassen multiply:", strassen(varC, varD))
|
||
|
||
# Example 3
|
||
int_print("Regular multiply:", varE %*% varF)
|
||
int_print("Strassen multiply:", strassen(varE, varF))
|
||
|
||
# Example 4
|
||
int_print("Regular multiply:", R %*% R)
|
||
int_print("Strassen multiply:", strassen(R, R))
|