This talk is about solving the linear system Cx = d when C is a Khatri-Rao (KR) product. This means that C is a block matrix where each block is a Kronecker product. Matrix-vector products that involve KR matrices can be computed very efficiently. Likewise, it is possible to solve a block triangular system efficiently if the matrix of coefficients is a KR product. Unfortunately, KR structure is lost when computing e.g., a Cholesky factorization. This prompts to develop an iterative framework for solving KR systems of equations. We develop the idea of an approximate KR factorization that can be used in conjunction with a Krylov method such as the conjugate gradient method. In particular, we show how to compute approximate Cholesky and LU factorizations with KR structure. The computation involves solving a nearest Kronecker product problem in each block position. Numerical results show that the approximate KR factorizations can significantly speed up the convergence of Krylov based methods when used as preconditioners.