snss_jd {SpatialBSS}

Spatial Non-Stationary Source Separation Joint Diagonalization

Description

snss_jd estimates the unmixing matrix assuming a spatial non-stationary source separation model implying non-constant covariance by jointly diagonalizing at least two covariance matrices computed for corresponding different sub-domains.

Usage

snss_jd(x, ...)

## Default S3 method:
snss_jd(x, coords, n_block, ordered = TRUE, ...)
## S3 method for class 'list'
snss_jd(x, coords, ordered = TRUE, ...)
## S3 method for class 'SpatialPointsDataFrame'
snss_jd(x, ...)
## S3 method for class 'sf'
snss_jd(x, ...)

Arguments

x

either a numeric matrix of dimension c(n, p) where the p columns correspond to the entries of the random field and the n rows are the observations, a list of length K defining the subdivision of the domain, an object of class sf or an object of class SpatialPointsDataFrame.

coords

a numeric matrix of dimension c(n,2) when x is a matrix where each row represents the sample location of a point in the spatial domain or a list of length K if x is a list which defines the subdivision of the domain. Not needed otherwise.

n_block

an integer defining the subdivision of the domain. See details.

ordered

logical. If TRUE the entries of the latent field are ordered by the sum of squared pseudo-eigenvalues of the diagonalized sub-domain covariance matrices. Default is TRUE.

...

further arguments for the fast real joint diagonalization algorithm that jointly diagonalizes the sub-domain covariance matrices. See details and frjd.

Details

This function assumes that the random field x is formed by

x(t) = A s(t) + b,

where A is the deterministic p \times p mixing matrix, b is the p-dimensional location vector, x is the observable p-variate random field given by the argument x, t are the spatial locations given by the argument coords and s is the latent p-variate random field assumed to consist of uncorrelated entries that have zero mean but non-constant variances. This function aims to recover s by

W(x(t) - \bar{x}),

where W is the p \times p unmixing matrix and \bar{x} is the sample mean. The function does this by splitting the given spatial domain into n_block^2 equally sized rectangular sub-domains and jointly diagonalizing the corresponding covariance matrices for all sub-domains.

Alternatively the domain subdivision can be defined by providing lists of length K for the arguments x and coords where the first list entries correspond to the values and coordinates of the first sub-domain and the second entries to the values and coordinates of the second sub-domain, etc..

snss_jd jointly diagonalizes the covariance matrices for each sub-domain with the function frjd. ... provides arguments for frjd, useful arguments might be:

• eps: tolerance for convergence.

• maxiter: maximum number of iterations.

Value

Similarly as sbss the function snss_jd returns a list of class 'snss' and 'sbss' with the following entries:

s

object of class(x) containing the estimated source random field.

coords

coordinates of the observations. Only given if x is a matrix or list.

w

estimated unmixing matrix.

w_inv

inverse of the estimated unmixing matrix.

d

matrix of stacked (jointly) diagonalized sub-domain covariance matrices with dimension c(n_block^2*p,p) or c(K*p,p) if x and coords are lists of length K.

x_mu

columnmeans of x.

cov_inv_sqrt

square root of the inverse sample covariance matrix of x.

References

Muehlmann, C., Bachoc, F. and Nordhausen, K. (2022), Blind Source Separation for Non-Stationary Random Fields, Spatial Statistics, 47, 100574, doi:10.1016/j.spasta.2021.100574.

See Also

sbss, sp, sf

Examples

# simulate coordinates
n <- 1000
coords <- runif(n * 2) * 20
dim(coords) <- c(n, 2)

# simulate random field
field_1 <- rnorm(n)
field_2 <- 2 * sin(pi / 20 * coords[, 1]) * rnorm(n)
field_3 <- rnorm(n) * (coords[, 1] < 10) + rnorm(n, 0, 3) * (coords[, 1] >= 10)

latent_field <- cbind(field_1, field_2, field_3)
mixing_matrix <- matrix(rnorm(9), 3, 3)
observed_field <- latent_field 

observed_field_sp <- sp::SpatialPointsDataFrame(coords = coords, 
                                                data = data.frame(observed_field))
sp::spplot(observed_field_sp, colorkey = TRUE, as.table = TRUE, cex = 1)

# apply snss_jd with 4 sub-domains 
res_4 <- snss_jd(observed_field, coords, n_block = 2)
JADE::MD(W.hat = coef(res_4), A = mixing_matrix)

# apply snss_jd with 9 sub-domains
res_9 <- snss_jd(observed_field, coords, n_block = 3)
JADE::MD(W.hat = coef(res_9), A = mixing_matrix)
cor(res_9$s, latent_field)

# print object
print(res_4)

# plot latent field
plot(res_4, colorkey = TRUE, as.table = TRUE, cex = 1)

# predict latent fields on grid
predict(res_4, colorkey = TRUE, as.table = TRUE, cex = 1)

# unmixing matrix
w_unmix <- coef(res_4)

# apply snss_jd with SpatialPointsDataFrame object 
res_4_sp <- snss_jd(observed_field_sp, n_block = 2)

# apply with list arguments
# first axis split by 5
# second axis split by 10
# results in 4 sub-domains
flag_x <- coords[, 1] < 5
flag_y <- coords[, 2] < 10
coords_list <- list(coords[flag_x & flag_y, ],
                    coords[!flag_x & flag_y, ],
                    coords[flag_x & !flag_y, ],
                    coords[!flag_x & !flag_y, ])
field_list <- list(observed_field[flag_x & flag_y, ],
                   observed_field[!flag_x & flag_y, ],
                   observed_field[flag_x & !flag_y, ],
                   observed_field[!flag_x & !flag_y, ])
plot(coords, col = 1)
points(coords_list[[2]], col = 2)
points(coords_list[[3]], col = 3)
points(coords_list[[4]], col = 4)

res_list <- snss_jd(x = field_list,
                    coords = coords_list)
plot(res_list, colorkey = TRUE, as.table = TRUE, cex = 1)
JADE::MD(W.hat = coef(res_list), A = mixing_matrix)