Schizophrenia is often associated with disrupted brain connectivity. However, identifying specific neuroimaging-based patterns pathognomonic for schizophrenia and related symptom severity remains a challenging open problem requiring large-scale datadriven analyses emphasizing not only statistical significance but also stability across multiple datasets, contexts and cohorts. Accurate prediction on previously unseen subjects, or generalization, is also essential for any useful biomarker of schizophrenia. In order to build a predictive model based on functional network feature patterns, we studied whole-brain fMRI functional networks, both at the voxel level and lower-resolution supervoxel level. Targeting Auditory Oddball task data on the FBIRN fMRI dataset (n = 95), we considered node-degree and link-weight network features and evaluated stability and generalization accuracy of statistically significant feature sets in discriminating patients vs. controls. We also applied sparse multivariate regression (elastic net) to whole-brain functional connectivity features, for the first time, to derive stable predictive features for symptom severity. Wholebrain link-weight features achieved 74% accuracy in identifying patients and were more stable than voxel-wise node-degrees. Linkweight features predicted severity of several negative and positive symptom scales, including inattentiveness and bizarre behavior. The most-significant, stable and discriminative functional connectivity changes involved increased correlations between thalamus and primary motor/primary sensory cortex, and between precuneus (BA7) and thalamus, putamen, and Brodmann areas BA9 and BA44. Precuneus, along with BA6 and primary sensory cortex, was also involved in predicting severity of several symptoms. Overall, the proposed multi-step methodology may help identify more reliable multivariate patterns allowing for accurate prediction of schizophrenia and its symptoms severity.