snss_sjd {SpatialBSS}

Spatial Non-Stationary Source Separation Spatial Joint Diagonalization

Description

snss_sjd estimates the unmixing matrix assuming a spatial non-stationary source separation model implying non-constant (spatial) covariance by jointly diagonalizing several covariance and/or spatial covariance matrices computed for a subdivision of the spatial domain into at least two sub-domains.

Usage

snss_sjd(x, ...)

## Default S3 method:
snss_sjd(x, coords, n_block, kernel_type = c('ring', 'ball', 'gauss'), 
     kernel_parameters, with_cov = TRUE, lcov = c('lcov', 'ldiff', 'lcov_norm'), 
     ordered = TRUE, ...)
## S3 method for class 'list'
snss_sjd(x, coords, kernel_type = c('ring', 'ball', 'gauss'), 
     kernel_parameters, with_cov = TRUE, lcov = c('lcov', 'ldiff', 'lcov_norm'), 
     ordered = TRUE, ...)
## S3 method for class 'SpatialPointsDataFrame'
snss_sjd(x, ...)
## S3 method for class 'sf'
snss_sjd(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

either be an integer defining the subdivision of the domain, 'x' or 'y'. See details.

kernel_type

a string indicating which kernel function to use. Either 'ring' (default), 'ball' or 'gauss'.

kernel_parameters

a numeric vector that gives the parameters for the kernel function. At least length of one for 'ball' and 'gauss' or two for 'ring' kernel, see details.

with_cov

logical. If TRUE not only spatial covariance matrices but also the sample covariances matrices for each sub-domain are considered in the joint diagonalization procedure. Default is TRUE.

lcov

a string indicating which type of local covariance matrix to use. Either 'lcov' (default), 'ldiff' or 'lcov_norm'. See sbss_asymp for details on the latter option.

ordered

logical. If TRUE the entries of the latent field are ordered by the sum of squared pseudo-eigenvalues of the diagonalized sub-domain (local) 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 (spatial) second order dependence. 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 spatial covariance matrices for all sub-domains. If the argument with_cov equals TRUE (default) then additionally also the sample covariance matrices for each sub-domain are included in the joint diagonalization procedure.

The arguments kernel_type, kernel_parameters and lcov determine which spatial kernel functions and which type of local covariance matrices are used for each sub-domain. The usage is equal to the function sbss.

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.. The argument n_block might be 'x' or 'y' indicating a split across the x or y coordinates similar as done by the function snss_sd.

snss_sjd 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 and/or local covariance matrices.

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_sjd with 4 sub-domains
# one ring kernel per sub-domain
# without covariances
res_4_ball <- snss_sjd(observed_field, coords, n_block = 2,
                  kernel_type = 'ball', kernel_parameters = c(0, 2), 
                  with_cov = TRUE)
JADE::MD(W.hat = coef(res_4_ball), A = mixing_matrix)

# apply snss_sjd with split across y
# one ring kernel per sub-domain
# without covariances
# should not work as field does not show spatial dependence
res_4_ring <- snss_sjd(observed_field, coords, n_block = 'y',
                       kernel_type = 'ring', kernel_parameters = c(0, 2), 
                       with_cov = FALSE)
JADE::MD(W.hat = coef(res_4_ring), A = mixing_matrix)

# print object
print(res_4_ball)

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

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

# unmixing matrix
w_unmix <- coef(res_4_ball)

# apply snss_jd with SpatialPointsDataFrame object 
res_4_ball_sp <- snss_sjd(observed_field_sp, n_block = 2,
                          kernel_type = 'ball', kernel_parameters = c(0, 2), 
                          with_cov = TRUE)

# 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_sjd(x = field_list,
                    coords = coords_list,
                    kernel_type = 'ring', kernel_parameters = c(0, 2))
plot(res_list, colorkey = TRUE, as.table = TRUE, cex = 1)
JADE::MD(W.hat = coef(res_list), A = mixing_matrix)