  The paper addresses joint sparsity selection in the regression coefficient
matrix and the error precision (inverse covariance) matrix for high-dimensional
multivariate regression models in the Bayesian paradigm. The selected sparsity
patterns are crucial to help understand the network of relationships between
the predictor and response variables, as well as the conditional relationships
among the latter. While Bayesian methods have the advantage of providing
natural uncertainty quantification through posterior inclusion probabilities
and credible intervals, current Bayesian approaches either restrict to specific
sub-classes of sparsity patterns and/or are not scalable to settings with
hundreds of responses and predictors. Bayesian approaches which only focus on
estimating the posterior mode are scalable, but do not generate samples from
the posterior distribution for uncertainty quantification. Using a bi-convex
regression based generalized likelihood and spike-and-slab priors, we develop
an algorithm called Joint Regression Network Selector (JRNS) for joint
regression and covariance selection which (a) can accommodate general sparsity
patterns, (b) provides posterior samples for uncertainty quantification, and
(c) is scalable and orders of magnitude faster than the state-of-the-art
Bayesian approaches providing uncertainty quantification. We demonstrate the
statistical and computational efficacy of the proposed approach on synthetic
data and through the analysis of selected cancer data sets. We also establish
high-dimensional posterior consistency for one of the developed algorithms.
