在scipy.optimize.minimizePython中为优化制定约束矩阵

如何解决在scipy.optimize.minimizePython中为优化制定约束矩阵

我正在尝试最小化以下功能：

设置是我想找到斜率和截距值，以使在轻度变化的光照条件下拍摄的许多图像上的像素值范围相等。因此，在这个问题中，矩阵D中使用了两种等式。

第一个方程是intercept + slope * pixel_val = measured_reflectance。因此，这里我们测量了一定波长范围内的真实表面反射率值。大约有6个这类方程。

在第二类方程式中，我们不知道表面反射率，而是尝试根据相机像素值和测得的表面反射率之间的线性关系进行估算。在这里，我们查看可以在两个重叠图像中看到的像素，并假设它们的值应该相同，这样 inter0+slope0*pixel1 == inter1+slope1*pixel2

因此，我将右侧移到左侧以将方程设置为零。这是针对在许多图像之间看到的大量像素完成的。

所以矩阵D代表一个系数矩阵，其中系数是图像中像素的值：

D = [
   [0,1,80],[1,45,64,0],14,10,...
   [0,52,50,0]
]

因此，每行表示一个像素，该像素在1张以上的图像中可见，1表示截距的系数，非1表示可以在图像中找到的特定像素的值。每列表示一个图片，因此对于上述D的第一第二行，它表明在图片1和图片3中看到了pixel1。

我正在尝试使用scipy.optimize.minimize函数执行此操作。我在弄清楚如何正确设置约束时遇到了麻烦。他们在论文中指出： “项Ax≤c表示集合中所有图像中所有像素的约束（作为估计系数的函数）。每个像素必须满足两个约束”

因此，我设置了两个约束函数（一个用于最小像素，一个用于最大像素）：

def cons_min(self,x):
    # constraint function for the minimum pixel values being above 0
    return np.array((x*self.C_min))
    
def cons_max(self,x):
    # constraint function for the maximum pixel values being below 1
    return np.array((x*self.C_max))-1

其中self.C_max和self.C_min都是形状为(68,136)的数组。由于约束被认为是估计系数的函数，因此我将最小和最大像素矩阵乘以x的系数。

这也是我的目标功能：

def f(self,x):
        # function to minimize
        y = np.dot(self.A,x) - self.B
        return (np.dot(y,y))*.5

self.A是形状为(31100,136)的稀疏矩阵/数组，而self.B是形状为(31100,)的数组。

最后，我使用以下方法来求解带约束的线性系统：

def solve(self):
    if self.A is None or self.B is None:
        m = "No equations found."
        raise AssertionError(m)

    initial = [1 if _ % 2 == 0 else 0 for _ in range(self.A.shape[1])]
        
    con1 = [{'type': 'ineq','fun': self.cons_min},{'type': 'ineq','fun': self.cons_max}
                ]
        
    res = optimize.minimize(self.f,initial,method='SLSQP',constraints=con1
                                )
        
    return res['x']

当前出现错误：

ValueError: all the input arrays must have same number of dimensions

我猜这是来自连接约束矩阵时出现的一个问题（c_ieq是我收集的不等式约束矩阵）。但是，这两个矩阵的维数应该相同，因此连接起来应该没有任何问题，所以我在这里出错的地方有点迷茫。

编辑： 这是我用来尝试实现此目的的完整代码。本质上，我有一个JSON文件，其中包含有关像素值，表面反射率以及图像之间具有相应像素值的匹配像素的信息。该代码设置了方程组，定义了目标函数和约束函数，然后具有最终求解方程的功能。

class SpectralRef():
    
    def __init__(self,matches,reflectances,mins,maxs):
        self.matches = matches
        self.reflectances = reflectances
        self.mins = mins
        self.maxs = maxs
        self.A = None
        self.B = None
        
        # collect image ids
        self.ids = list(matches.keys())
        for panel in reflectances:
            for image in reflectances[panel]:
                if image not in self.ids:
                    self.ids.append(image)
        self.valid = set()

    def equations(self):
        """Provides the equations to fit the model.
    
        Parameters
        ----------
        image_means : dict(Number)
            Mean values of all images.
        reflectances : optional,dict(dict({"values": list(Number),"reflectance": list(Number)}))
            Digital number and reflectance value pair for some of the images.
        min_matches : optional,positive int
            Minimum matches to use an image for creating the equation system.
    
        Returns
        -------
        A : np.ndarray(Number,shape=(n,k))
            Matrix of equations.
        B : np.ndarray(Number,shape=(n))
            Vector of solutions.
    
        Examples
        --------
        
        Create a matches dictionary for a single band.
        matches = {'TCR00109_0': 
                    {'TCR00121_0': {
                        'values_r': [21376,17856,24512,16256]
                        'values_n': [13440,8064,17728,8576]
                        'distance_center_r': [[466.59501922,180.202608],[437.16801672,412.052254],[418.18338466,195.249375],[399.39276685,125.824875]
                                             ]
                        'distance_center_n' : [[466.59501922,125.824875]
                                             ]
                                    }
                    }
                }

    
    
        """
        n = len(self.ids)
    
        # validate input
        if not isinstance(self.matches,dict):
            m = "'matches' needs to be a dict,got %s"
            raise TypeError(m % type(self.matches))
        for id_r in self.matches:
            if not isinstance(self.matches[id_r],dict):
                m = "all items of 'matches' needs to be a dict,got %s for %s"
                raise TypeError(m % (type(self.matches[id_r]),id_r))
            for id_n in self.matches[id_r]:
                if id_n not in self.ids:
                    m = "Unknown image id %s"
                    raise ValueError(m % id_n)
                    
                #check structure of sub-dictionaries
                pairs = self.matches[id_r][id_n]
    
                # convert to dict,if necessary
                if isinstance(pairs,np.recarray):
                    pairs = {key: pairs[key] for key in pairs.dtype.names}
    
                if not isinstance(pairs,dict):
                    m = "'pairs' needs to be a dictionary,got %s"
                    raise TypeError(m % type(pairs))
    
                if "values_r" not in pairs:
                    m = "Dictionary key 'values_r' not found in %s %s pair"
                    raise KeyError(m % (id_r,id_n))
                if "values_n" not in pairs:
                    m = "Dictionary key 'values_n' not found in %s %s pair"
                    raise KeyError(m % (id_r,id_n))
                pairs["values_r"] = np.array(pairs["values_r"])
                pairs["values_n"] = np.array(pairs["values_n"])
                if not len(pairs["values_r"]) == len(pairs["values_n"]):
                    m = "length of 'values_r' and 'values_n' differ for %s-%s"
                    raise ValueError(m % (id_r,id_n))
    
                if "weights" in pairs:
                    pairs["weights"] = np.array(pairs["weights"])
                    if not len(pairs["weights"]) == len(pairs["values_r"]):
                        m = "length of 'weights' and 'values_n' differ for %s-%s"
                        raise ValueError(m % (id_r,id_n))
    
        # Define equation sysem Ax = B
        A = []
        B = []
    
        print('add reflectance equations')
        
        # add reflectance equations
        for panel in self.reflectances:
            # Only take the most nadir images from each panel as the reflectance
            for image in self.reflectances[panel]:
                img_id = self.ids.index(image)
                a = np.zeros((1,2 * n),dtype=np.float64)
                d = self.reflectances[panel][image]
                a[:,2*img_id] = d['values']
                a[:,2*img_id + 1] = 1
                A.append(a)
        
                b = np.zeros(1,dtype=np.float64)
                b[:] = d['reflectances']
                B.append(b)

        # Maybe add in equations where each reference panel is seen in multiple
        # images. Then we could make these equations similar to the tiepoints
        # where the reflectance of the panel in image 1 should equal the reflect
        # of the same panel in image 2. Try this if the first approach is not 
        # very good
    
        # add pairs
        print('add tiepoint pair equations')
        for id_r in self.matches:
            # Skip if not in the list of ids
            if id_r not in self.ids:
                continue
            for id_n in self.matches[id_r]:
                
                if id_n not in self.ids:
                    continue
                if id_r == id_n:
                    continue
                # find indices
                idx_r = self.ids.index(id_r)
                idx_n = self.ids.index(id_n)
    
                # get all matches
                d = self.matches[id_r][id_n]
                
                # Skip if less than 3 tiepoints
                if len(d['values_r']) < 20:
                    continue
                
                for vr,vn in zip(d["values_r"],d["values_n"]):
                    a = np.zeros((1,dtype=np.float64)
                    # Try swapping the (-),but I doubt the order of this matters
                    a[0,2*idx_r] = vr
                    a[0,2*idx_r+1] = 1 
                    a[:,2*idx_n] = -vn
                    a[:,2*idx_n+1] = -1
                    b = np.zeros(1,dtype=np.float64)
        
                    A.append(a)
                    B.append(b)
    
        # add image min/max values
        A = np.vstack(A)
        B = np.concatenate(B)
        self.A = A
        self.B = B
        
    def constraints(self):
        C_min = []
        C_max = []
        n = len(self.ids)
        
        for id_r in self.ids:
            a = np.zeros((1,dtype=np.float64)
            idx_r = self.ids.index(id_r)
            value = self.maxs[id_r]
            a[0,2*idx_r] = -value
            a[0,2*idx_r + 1] = -1    
            C_max.append(a)

        for id_r in self.ids:
            a = np.zeros((1,dtype=np.float64)
            idx_r = self.ids.index(id_r)
            value = self.mins[id_r]
            a[0,2*idx_r] = value
            a[0,2*idx_r + 1] = 1
            C_min.append(a)
            
        self.C_min = np.vstack(C_min)
        self.C_max = np.vstack(C_max)
        
    def f(self,y))*.5
        
    def cons_min(self,x):
        # constraint function for the minimum pixel values
        return np.array((x*self.C_min))
    
    def cons_max(self,x):
        # constraint function for the maximum pixel values
        return np.array((x*self.C_max))-1
    
    def solve(self):
        if self.A is None or self.B is None:
            m = "No equations found."
            raise AssertionError(m)

        

        initial = [1 if _ % 2 == 0 else 0 for _ in range(self.A.shape[1])]
        
        con1 = [{'type': 'ineq','fun': self.cons_max}
                ]
        
        res = optimize.minimize(self.f,constraints=con1
                                )
        
        return res['x']