Suppose I have the following quadratic program
$$\begin{array}{ll} \underset{v_1, v_2, \dots, v_m \in \mathbb{R}^n}{\text{minimize}} & \displaystyle\sum_i v_i^T a_i + \displaystyle\sum_{i,j} \lambda_{i,j} \, v_i^T A v_j\\ \text{subject to} & \displaystyle\sum_i v_i^T b_{k,i} = c_k, \quad \forall k\\ & v_1, v_2, \dots, v_m \geq \vec{0}\end{array}$$
where $a_i, b_{k,i} \in \mathbb{R}^n$, $\lambda_{i,j}, c_k \in \mathbb{R}$, and $A \in \mathbb{R}^{n \times n}$.
How to solve a quadratic program like this? Is there any way to re-write in canonical QP form? If not, is there any Python package that can solve quadratic program like this?
(edit) Is $\sum_{i,j}\lambda_{i,j} \: v_i^TAv_j$ just
$$\begin{bmatrix} v_1^T &v_2^T & \dotsb & v_m^T\end{bmatrix} \begin{bmatrix} \lambda_{1,1}A&\lambda_{1,2}A& \dotsb & \lambda_{1,m}A\\\lambda_{2,1}A&\lambda_{2,2}A& \dotsb & \lambda_{2,m}A\\ \vdots&\vdots& \ddots & \lambda_{m-1,m}A \\\lambda_{m,1}A&\lambda_{m,2}A & \dotsb &\lambda_{m,m}A\end{bmatrix}\begin{bmatrix} v_1^T \\v_2^T \\ \vdots \\ v_m^T\end{bmatrix} $$
So it's equivalent to the one variable situation?
Let's set
$$ \begin{align} \bf{v} &= (v_1, \ldots, v_m) = (v_{11}, \ldots, v_{1n}, \ldots, v_{m1}, \ldots, v_{mn}) \in \mathbb{R}^{n \cdot m} \\ % \bf{a} &= (a_1, \ldots, a_m) = (a_{11}, \ldots, a_{1n}, \ldots, a_{m1}, \ldots, a_{mn}) \in \mathbb{R}^{n \cdot m} \\ % \bf{B} &= \begin{pmatrix} % b_{111} & \ldots & b_{11n} & \ldots & b_{1m1} & \ldots & b_{1mn} \\ \vdots & & & & & & \vdots \\ b_{p11} & \ldots & b_{p1n} & \ldots & b_{pm1} & \ldots & b_{pmn}\\ \end{pmatrix} \in \mathbb{R}^{p \times m\cdot n} \\ % c &= (c_1, \ldots, c_p) \in \mathbb{R}^{p} \\ \Lambda &= (\lambda_{i,j})_{i = 1, \ldots, n, j = 1, \ldots, n} \\ {\bf D} &= \text{diag}{(\Lambda \cdot A, \ldots, \Lambda \cdot A}) \in \mathbb{R}^{m \cdot n \times m \cdot n} \end{align} $$
Then, we have
$$ \begin{align} \sum_{i=1}^{m} v_i^{\top} b_{k,i} &= c_k \quad \forall k = 1,\ldots, p \iff {\bf B v} = c. \\ % \sum_{i=1}^{m} v_i^{\top} a_i &= \bf{a}^{\top} v \\ \sum_{i=1}^{m} \sum_{i=1}^{m} \lambda_{i,j} v_i^{\top} A v_j &= {\bf{v^{\top}}} {\bf D} {\bf v} \\ \end{align} $$
and consequently
$$ \begin{align} &\min_{{\bf v}} &\quad {\bf{v^{\top}}} {\bf D} {\bf v} + {\bf{a}}^{\top} {\bf v} \\ &\text{s.t.} &\quad {\bf B v} &= c \\ && {\bf v} &\geq 0 \end{align} $$
However, as already mentioned in the comments, non-convex QPs are quite hard to solve towards global optimality. Therefore, one usually settles for a local minimizer. Here, you can use the Ipopt solver (cyipopt is a python interface for the Ipopt solver). In case you really need a global minimum, you can use the commercial solver Gurobi. They offer free academic licenses.