Multispectral　bioluminescence　tomography　is　becoming　a　promising　tool　because　it　can　resolve　the　biodistibution　of　bioluminescent　reporters　associated　with　cellular　and　subcellular　function　through　several　millimeters　with　to　centimeters　of　tissues　in　vivo.　Generally,　to　recover　the　bioluminescent　sources,　the　source　reconstruction　problem　is　formulated　as　a　nonlinear　least-squares-type　bounds　constrained　optimization　problem.　However,　bioluminescence　tomography　(BLT)　is　an　ill-posed　problem.　For　the　sake　of　stability　and　uniqueness　of　BLT,　many　algorithms　have　been　proposed　to　regularize　the　problem,　such　as　L-2　norm　and　L-1　norm.　Here,　we　proposed　a　new　regularization　method　with　Huber　function　to　regularize　BLT　problem　to　obtain　robustness　like　L-1　and　rapid　convergence　of　L-2.　Furthermore,　the　computational　burden　is　largely　increased　with　the　use　of　spectral　data.　Therefore,　there　is　a　critical　need　to　develop　a　fast　reconstruction　algorithm　for　solving　multispectral　bioluminescence　tomography.　In　the　paper,　a　limited　memory　quasi-Newton　algorithm　for　solving　the　large-scale　optimization　problem　is　proposed　to　fast　localize　the　bioluminescent　source.　In　the　numerical　simulation,　a　heterogeneous　phantom　was　used　to　evaluate　the　performance　of　the　proposed　algorithm　with　the　Monte　Carlo　based　synthetic　data.　Additionally,　the　real　mouse　experiments　were　conducted　to　further　evaluate　the　proposed　algorithm.　The　results　demonstrate　the　potential　and　merits　of　the　proposed　algorithm.