This is a common problem in covariance matrix estimation, with several possible solutions. One of the simplest involves two steps:
(1) You compute each element of the covariance matrix on a 'best efforts' basis, meaning you take the covariance of the two time series involved after REMOVING any data pairs having a N/A value. (Note that this means each element of the matrix will be based on a different number of observations, which means the resulting matrix is not a standard covariance matrix, it may not be positive definite for example). I assume that for any two time series there are at least a few common observations, otherwise it is an ill-posed problem as Matthew Gunn pointed out.
(2) You "massage" the resulting matrix to make it positive definite (and thus acceptable for use as a covariance matrix) using the routine nearPD which is available in R link
[Even after all this work, a large covariance matrix will be very 'noisy' and of poor quality. You should consider further steps such as 'shrinkage' link before you use the results to find an optimum portfolio].