Accurately predicting response properties of molecules such as the dynamic polarizability and hyperpolarizability using quantum mechanics has been a long-standing challenge with widespread applications in material and drug design. Classical simulation techniques in quantum chemistry are hampered by the exponential growth of the many-electron Hilbert space as the system size increases... In this work, we propose an algorithm for computing linear and nonlinear molecular response properties on quantum computers, by first reformulating the target property into a symmetric expression more suitable for quantum computation via introducing a set of auxiliary quantum states, and then determining these auxiliary states via solving the corresponding linear systems of equations on quantum computers. On one hand, we prove that using the quantum linear system algorithm [Harrow et al., Phys. Rev. Lett. 103, 150502 (2009)] as a subroutine the proposed algorithm scales only polynomially in the system size instead of the dimension of the exponentially large Hilbert space, and hence achieves an exponential speedup over existing classical algorithms. On the other hand, we introduce a variational hybrid quantum-classical variant of the proposed algorithm, which is more practical for near-term quantum devices. (read more)