Traditional task-based fMRI activation detection methods, e.g., the widely used general linear model (GLM), assume that the brain's hemodynamic responses follow the block-based or event-related stimulus paradigm. Typically, these activation detections are performed voxel-wise independently, and then are usually followed by statistical corrections. Despite remarkable successes and wide adoption of these methods, it remains largely unknown how functional brain regions interact with each other within specific networks during task performance blocks and in the baseline. In this paper, we present a novel algorithmic pipeline to statistically infer and sparsely represent higher-order functional interaction patterns within the working memory network during task performance and in the baseline. Specifically, a collection of higher-order interactions are inferred via the greedy equivalence search (GES) algorithm for both task and baseline blocks. In the next stage, an effective online dictionary learning algorithm is utilized for sparse representation of the inferred higher-order interaction patterns. Application of this framework on a working memory task-based fMRI data reveals interesting and meaningful distributions of the learned sparse dictionary atoms in task and baseline blocks. In comparison with traditional voxel-wise activation detection and recent pair-wise functional connectivity analysis, our framework offers a new methodology for representation and exploration of higher-order functional activities in the brain.