D.matrix and Duplication.matrix function in 'R' package 'matrixcalc'

Win_odd Dhamnekar

Win_odd Dhamnekar

Junior Member
Joined
Aug 14, 2018
Messages
212
In 'R' while using D.matrix(n) or duplication.matrix (n) function,where n=3, we get following answer.
> D.matrix(3)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 0 0 0 0 0
[2,] 0 1 0 0 0 0
[3,] 0 0 1 0 0 0
[4,] 0 1 0 0 0 0
[5,] 0 0 0 1 0 0
[6,] 0 0 0 0 1 0
[7,] 0 0 1 0 0 0
[8,] 0 0 0 0 1 0
[9,] 0 0 0 0 0 1
>
Where

> u.vectors(3)
$k
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 2 4 0
[3,] 3 5 6

$I
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 0 0 0 0 0
[2,] 0 1 0 0 0 0
[3,] 0 0 1 0 0 0
[4,] 0 0 0 1 0 0
[5,] 0 0 0 0 1 0
[6,] 0 0 0 0 0 1
Source programme(code) for D.matrix is below.
> getAnywhere(D.matrix)
A single object matching ‘D.matrix’ was found
It was found in the following places
package:matrixcalc
namespace:matrixcalc
with value

function (n)
{
if (missing(n))
stop("argument n is missing")
if (!is.numeric(n))
stop("argument n is not numeric")
if (n != trunc(n))
stop("argument n is not an integer")
if (n < 2)
stop("argument n is less than 2")
p <- n * (n + 1)/2
nsq <- n * n
Dt <- matrix(0, nrow = p, ncol = nsq)
T <- T.matrices(n)
u <- u.vectors(n)
k <- u$k
I <- u$I
for (j in 1:n) {
for (i in j:n) {
Dt <- Dt + I[, k[i, j]] %*% t(vec(T[][[j]]))
}
}
return(t(Dt))
}
<bytecode: 0x000002926d7892e0>
<environment: namespace:matrixcalc>
> T.matrices(3)
[[1]]
[[1]][[1]]
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 0 0
[3,] 0 0 0

[[1]][[2]]
[,1] [,2] [,3]
[1,] 0 1 0
[2,] 1 0 0
[3,] 0 0 0

[[1]][[3]]
[,1] [,2] [,3]
[1,] 0 0 1
[2,] 0 0 0
[3,] 1 0 0


[[2]]
[[2]][[1]]
[,1] [,2] [,3]
[1,] 0 1 0
[2,] 1 0 0
[3,] 0 0 0

[[2]][[2]]
[,1] [,2] [,3]
[1,] 0 0 0
[2,] 0 1 0
[3,] 0 0 0

[[2]][[3]]
[,1] [,2] [,3]
[1,] 0 0 0
[2,] 0 0 1
[3,] 0 1 0


[[3]]
[[3]][[1]]
[,1] [,2] [,3]
[1,] 0 0 1
[2,] 0 0 0
[3,] 1 0 0

[[3]][[2]]
[,1] [,2] [,3]
[1,] 0 0 0
[2,] 0 0 1
[3,] 0 1 0

[[3]][[3]]
[,1] [,2] [,3]
[1,] 0 0 0
[2,] 0 0 0
[3,] 0 0 1

Now I don't understand how did this 'matrixcalc' package compute t ( vec ( T [ [ i ]] [ [ j ] ] ) ) and multiply itself with Dt + I [ , k[ i , j ] ]?
Would any member of this forum explain?

I shall ask this question to author of this package also.
 
Please the read the nested loop as follows :
for (j in 1:n) {
for (i in j:n) {
Dt <- Dt + I [ , k[ i , j ] ] %*% t( vec( T [ [ i ] ] [ [ j ] ] ) )
}
}
 
I got 9 x 6 t (vec ( T [ [ i] ] [ [ j ] ] ) ) matrix by doing following computations in 'R'.

> # Create a 6x6 identity matrix
> identity_matrix <- diag(6)
>
> # Print the identity matrix
> print(identity_matrix)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 0 0 0 0 0
[2,] 0 1 0 0 0 0
[3,] 0 0 1 0 0 0
[4,] 0 0 0 1 0 0
[5,] 0 0 0 0 1 0
[6,] 0 0 0 0 0 1
> matrix.inverse(identity_matrix)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 0 0 0 0 0
[2,] 0 1 0 0 0 0
[3,] 0 0 1 0 0 0
[4,] 0 0 0 1 0 0
[5,] 0 0 0 0 1 0
[6,] 0 0 0 0 0 1
> identity_matrix %*% d
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 1 0 0 0 0 0 0 0 0
[2,] 0 1 0 1 0 0 0 0 0
[3,] 0 0 1 0 0 0 1 0 0
[4,] 0 0 0 0 1 0 0 0 0
[5,] 0 0 0 0 0 1 0 1 0
[6,] 0 0 0 0 0 0 0 0 1
> t(identity_matrix %*% d)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 0 0 0 0 0
[2,] 0 1 0 0 0 0
[3,] 0 0 1 0 0 0
[4,] 0 1 0 0 0 0
[5,] 0 0 0 1 0 0
[6,] 0 0 0 0 1 0
[7,] 0 0 1 0 0 0
[8,] 0 0 0 0 1 0
[9,] 0 0 0 0 0 1
> D.matrix(3) == t(identity_matrix %*% d)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] TRUE TRUE TRUE TRUE TRUE TRUE
[2,] TRUE TRUE TRUE TRUE TRUE TRUE
[3,] TRUE TRUE TRUE TRUE TRUE TRUE
[4,] TRUE TRUE TRUE TRUE TRUE TRUE
[5,] TRUE TRUE TRUE TRUE TRUE TRUE
[6,] TRUE TRUE TRUE TRUE TRUE TRUE
[7,] TRUE TRUE TRUE TRUE TRUE TRUE
[8,] TRUE TRUE TRUE TRUE TRUE TRUE
[9,] TRUE TRUE TRUE TRUE TRUE TRUE

But now, if I want to prepare t(vec( T [ [ i ] ] [ [ j ] ] ) ) , how can I prepare it in 'R'?
 
Please note that first I would need to prepare manually T[ [i ] [ [ j ] ] matrix. Then I need to vectorize it by vec ( ) command. Then, transposing the resultant matrix, I get t ( vec ( T [ [ i ] ] [ [ j ] ] ) ). The question is how to make T [ [ i ] [ [ j ] ] in 'R' ?
 
I suspect that most, myself included, helpers here don't know R. You would probably have better luck in getting help if your question could be reduced to plain math.
 
You must log in or register to reply here.
Top