Fractals and Kronecker product

[This article was first published on Freakonometrics - Tag - R-english, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

A few years ago, I went to listen to Roger Nelsen who was giving a talk about copulas with fractal support. Roger is amazing when he gives a talk (I am also a huge fan of his books, and articles), and I really wanted to play with that concept (that he did publish later on, with Gregory Fredricks and José Antonio Rodriguez-Lallena). I did mention that idea in a paper, writen with Alessandro Juri, just to mention some cases where deriving fixed point theorems is not that simple (since the limit may not exist).

The idea in the initial article was to start with something quite simple, a the so-called transformation matrix, e.g.

http://latex.codecogs.com/gif.latex?T=\frac{1}{8}\left(\begin{matrix}1&%200%20&%201%20\\%200%20&%204%20&%200%20\\%201%20&%200&1\end{matrix}\right)

Here, in all areas with mass, we spread it uniformly (say), i.e. the support of http://latex.codecogs.com/gif.latex?T(C^\perp) is the one below, i.e. http://latex.codecogs.com/gif.latex?1/8th of the mass is located in each corner, and http://latex.codecogs.com/gif.latex?1/2 is in the center. So if we spread the mass to have a copula (with uniform margin,)we have to consider squares on intervals http://latex.codecogs.com/gif.latex?[0,1/4], http://latex.codecogs.com/gif.latex?[1/4,3/4] and http://latex.codecogs.com/gif.latex?[3/4,1],

Then the idea, then, is to consider http://latex.codecogs.com/gif.latex?T^2=\otimes^2T, where  http://latex.codecogs.com/gif.latex?\otimes^2T is the tensor product (also called Kronecker product) of http://latex.codecogs.com/gif.latex?T with itself. Here, the support of http://latex.codecogs.com/gif.latex?T^2(C^\perp) is

Then, consider http://latex.codecogs.com/gif.latex?T^3=\otimes^3T, where http://latex.codecogs.com/gif.latex?\otimes^3T is the tensor product of http://latex.codecogs.com/gif.latex?T with itself, three times. And the support of http://latex.codecogs.com/gif.latex?T^3(C^\perp) is

Etc. Here, it is computationally extremely simple to do it, using this Kronecker product. Recall that if http://latex.codecogs.com/gif.latex?%20%20%20%20%20\mathbf{A}=(a_{i,j}), then

http://latex.codecogs.com/gif.latex?%20%20%20%20%20\mathbf{A}\otimes\mathbf{B}%20=%20\begin{pmatrix}%20a_{11}%20\mathbf{B}%20&%20\cdots%20&%20a_{1n}\mathbf{B}%20\\%20\vdots%20&%20\ddots%20&%20\vdots%20\\%20a_{m1}%20\mathbf{B}%20&%20\cdots%20&%20a_{mn}%20\mathbf{B}%20\end{pmatrix}

So, we need a transformation matrix: consider the following http://latex.codecogs.com/gif.latex?4\times4 matrix,
> k=4
> M=matrix(c(1,0,0,1,
+            0,1,1,0,
+            0,1,1,0,
+            1,0,0,1),k,k)
> M
[,1] [,2] [,3] [,4]
[1,]    1    0    0    1
[2,]    0    1    1    0
[3,]    0    1    1    0
[4,]    1    0    0    1
Once we have it, we just consider the Kronecker product of this matrix with itself, which yields a http://latex.codecogs.com/gif.latex?4^2\times4^2 matrix,
> N=kronecker(M,M)
> N[,1:4]
[,1]  [,2] [,3] [,4]
[1,]     1    0    0    1
[2,]     0    1    1    0
[3,]     0    1    1    0
[4,]     1    0    0    1
[5,]     0    0    0    0
[6,]     0    0    0    0
[7,]     0    0    0    0
[8,]     0    0    0    0
[9,]     0    0    0    0
[10,]    0    0    0    0
[11,]    0    0    0    0
[12,]    0    0    0    0
[13,]    1    0    0    1
[14,]    0    1    1    0
[15,]    0    1    1    0
[16,]    1    0    0    1
And then, we continue,
> for(s in 1:3){N=kronecker(N,M)}
After only a couple of loops, we have a http://latex.codecogs.com/gif.latex?4^5\times4^5 matrix. And we can plot it simply to visualize the support,
> image(N,col=c("white","blue"))
As we zoom in, we can visualize this fractal property,

To leave a comment for the author, please follow the link and comment on their blog: Freakonometrics - Tag - R-english.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)