Is it accurate compared to bivariate JS-divergence?
Yes it is.
Would a matrix of bivariate JS-divergences feasible (like the correlation matrix), and what would its diagonal consist of?
It is feasible, I made one once. the diagonal would be equal to 0 because divergence is null between pairs of the same distribution.
Anyway, if you have a matrix of Jensen-Shannon divergences, you are probably doing clustering or some distance-based analysis anyway. In that case, you may prefer to take the square root of the divergences. What you get then is called Jensen-Shannon distance.
Besides these, if t-distributed Stochastic Neighbor Embedding (t-SNE) in sklearn.manifold.tsne can be used to form a representation of multivariate KL-divergence for minimization (multivariate rather than one pair at a time), can the same or another algorithm do the same for a multivariate version of Jensen-Shannon divergence?
For what I know, there is no multivariate KL divergence involved in t-SNE algorithm. KL divergence is minimized between the original distribution od the distances and the one following dimensionality reduction. In general, it makes more sense to use KL divergence when you have some given distribution and you want to approximate it with a simpler one, and use JS divergence instead when you have two distributions that are both taken from data and therefore given. In conclusion, you may find yourself minimizing KL divergence, but I can't find any reasonable setting where to minimize JS divergence. This nice answer elaborates in this issue, with some good insight.