189 8069 5689
验证码（CAPTCHA）全称为全自动区分计算机和人类的公开图灵测试（Completely Automated Public Turing test to tell Computersand Humans Apart）。从其全称可以看出，验证码用于测试用户是真实的人类还是计算机机器人。
在红山等地区，都构建了全面的区域性战略布局，加强发展的系统性、市场前瞻性、产品创新能力，以专注、极致的服务理念，为客户提供网站建设、成都做网站 网站设计制作按需规划网站,公司网站建设,企业网站建设,品牌网站建设,全网营销推广,外贸网站建设,红山网站建设费用合理。
1.获得验证码图片
每次加载注册网页都会显示不同的验证验图像，为了了解表单需要哪些参数，我们可以复用上一章编写的parse_form()函数。
import cookielib,urllib2,pprint import form REGISTER_URL = '' cj=cookielib.CookieJar() opener=urllib2.build_opener(urllib2.HTTPCookieProcessor(cj)) html=opener.open(REGISTER_URL).read() form=form.parse_form(html) pprint.pprint(form)
{'_formkey': 'a67cbc84-f291-4ecd-9c2c-93937faca2e2', '_formname': 'register', '_next': '/places/default/index', 'email': '', 'first_name': '', 'last_name': '', 'password': '', 'password_two': '', 'recaptcha_response_field': None} 123456789101112131415161718
上面recaptcha_response_field是存储验证码的值，其值可以用Pillow从验证码图像获取出来。先安装pip install Pillow，其它安装Pillow的方法可以参考 。Pillow提价了一个便捷的Image类，其中包含了很多用于处理验证码图像的高级方法。下面的函数使用注册页的HTML作为输入参数，返回包含验证码图像的Image对象。
import lxml.html from io import BytesIO from PIL import Image tree=lxml.html.fromstring(html) print tree
Element html at 0x7f8b006ba890 img_data_all=tree.cssselect('div#recaptcha img')[0].get('src') print img_data_all
data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAQAAAABgCAIAAAB9kzvfAACAtklEQVR4nO29Z5gcZ5ku3F2dc865
...
rkJggg== img_data=img_data_all.partition(',')[2] print img_data
iVBORw0KGgoAAAANSUhEUgAAAQAAAABgCAIAAAB9kzvfAACAtklEQVR4nO29Z5gcZ5ku3F2dc865
...
rkJggg== binary_img_data=img_data.decode('base64') file_like=BytesIO(binary_img_data) print file_like
_io.BytesIO object at 0x7f8aff6736b0 img=Image.open(file_like) print img
PIL.PngImagePlugin.PngImageFile image mode=RGB size=256x96 at 0x7F8AFF5FAC90 12345678910111213141516171819202122232425
在本例中，这是一张进行了Base64编码的PNG图像，这种格式会使用ASCII编码表示二进制数据。我们可以通过在第一个逗号处分割的方法移除该前缀。然后，使用Base64解码图像数据，回到最初的二进制格式。要想加载图像，PIL需要一个类似文件的接口，所以在传给Image类之前，我们以使用了BytesIO对这个二进制数据进行了封装。
完整代码:
# -*- coding: utf-8 -*-form.pyimport urllibimport urllib2import cookielibfrom io import BytesIOimport lxml.htmlfrom PIL import Image
REGISTER_URL = ''#REGISTER_URL = ''def extract_image(html):
tree = lxml.html.fromstring(html)
img_data = tree.cssselect('div#recaptcha img')[0].get('src') # remove data:image/png;base64, header
img_data = img_data.partition(',')[-1] #open('test_.png', 'wb').write(data.decode('base64'))
binary_img_data = img_data.decode('base64')
file_like = BytesIO(binary_img_data)
img = Image.open(file_like) #img.save('test.png')
return imgdef parse_form(html):
"""extract all input properties from the form
"""
tree = lxml.html.fromstring(html)
data = {} for e in tree.cssselect('form input'): if e.get('name'):
data[e.get('name')] = e.get('value') return datadef register(first_name, last_name, email, password, captcha_fn):
cj = cookielib.CookieJar()
opener = urllib2.build_opener(urllib2.HTTPCookieProcessor(cj))
html = opener.open(REGISTER_URL).read()
form = parse_form(html)
form['first_name'] = first_name
form['last_name'] = last_name
form['email'] = email
form['password'] = form['password_two'] = password
img = extract_image(html)#
captcha = captcha_fn(img)#
form['recaptcha_response_field'] = captcha
encoded_data = urllib.urlencode(form)
request = urllib2.Request(REGISTER_URL, encoded_data)
response = opener.open(request)
success = '/user/register' not in response.geturl() #success = '/places/default/user/register' not in response.geturl()
return success12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152
2.光学字符识别验证码
光学字符识别（Optical Character Recognition, OCR）用于图像中抽取文本。本节中，我们将使用开源的Tesseract OCR引擎，该引擎最初由惠普公司开发的，目前由Google主导。Tesseract的安装说明可以从 获取。然后可以使用pip安装其Python封装版本pytesseractpip install pytesseract。
下面我们用光学字符识别图像验证码：
import pytesseract import form img=form.extract_image(html) pytesseract.image_to_string(img)'' 123456
如果直接把验证码原始图像传给pytesseract，一般不能解析出来。这是因为Tesseract是抽取更加典型的文本，比如背景统一的书页。下面我们进行去除背景噪音，只保留文本部分。验证码文本一般都是黑色的，背景则会更加明亮，所以我们可以通过检查是否为黑色将文本分离出来，该处理过程又被称为阈值化。
img.save('2captcha_1original.png') gray=img.convert('L') gray.save('2captcha_2gray.png') bw=gray.point(lambda x:0 if x1 else 255,'1') bw.save('2captcha_3thresholded.png') 1234567
这里只有阈值小于1的像素（全黑）都会保留下来，分别得到三张图像：原始验证码图像、转换后的灰度图和阈值化处理后的黑白图像。最后我们将阈值化处理后黑白图像再进行Tesseract处理，验证码中的文字已经被成功抽取出来了。
pytesseract.image_to_string(bw)'language' import Image,pytesseract img=Image.open('2captcha_3thresholded.png') pytesseract.image_to_string(img)'language' 123456789
我们通过示例样本测试，100张验证码能正确识别出90张。
import ocr ocr.test_samples()
Accuracy: 90/100 1234
下面是注册账号完整代码：
# -*- coding: utf-8 -*-import csvimport stringfrom PIL import Imageimport pytesseractfrom form import registerdef main():
print register('Wu1', 'Being1', 'Wu_Being001@qq.com', 'example', ocr)def ocr(img):
# threshold the image to ignore background and keep text
gray = img.convert('L') #gray.save('captcha_greyscale.png')
bw = gray.point(lambda x: 0 if x 1 else 255, '1') #bw.save('captcha_threshold.png')
word = pytesseract.image_to_string(bw)
ascii_word = ''.join(c for c in word if c in string.letters).lower() return ascii_wordif __name__ == '__main__':
main()1234567891011121314151617181920212223
我们可以进一步改善OCR性能：
- 实验不同阈值
- 腐蚀阈值文本，突出字符形状
- 调整图像大小
- 根据验证码字体训练ORC工具
- 限制结果为字典单词
一般模型既有不等式约束，也有等式约束；既有非负的约束决策变量，也有整个实数域上的自由决策变量。
标准模型
引入冗余的决策变量，使得不等式约束转化为等式约束。这里的每个决策变量都具有非负性。
在这里插入图片描述
把上述模型用矩阵表示就是
m i n ( o r m a x ) C T X s . t A X = b ⃗ X ≥ 0 min(or\ max) C^TX\\ s.t \ AX=\vec{b}\\ \ X \geq 0
min(or max)C
T
X
s.t AX=
b
X≥0
线性规划问题的基本假设
系数矩阵A的行向量线性无关。
如果线性相关有2种可能，要么是增广矩阵的该行也线性相关，则该行约束冗余，可以删去。要么增广矩阵的该行线性无关，则方程无解，优化问题不存在。
系数矩阵A的行数小于列数
如果行数m大于列数n，则行向量是m个n维向量，不可能线性无关。吐过行数等于列数，且行向量线性无关，则约束条件确定了唯一解，无需优化。
一般模型与标准模型的转化
主要方式是增加决策变量。有两种情况需要增加
不等式变等式，每个不等式增加一个决策变量。
1个自由决策变量转化为2个约束的决策变量。
在这里插入图片描述
线性规划问题解的可能情况
唯一最优解
没有有限的最优目标函数
没有可行解
无穷多的最优解（一维问题中不会出现）
凸集
Def. 凸集：该集合中任意两个元素的凸组合仍然属于该集合。
在这里插入图片描述
注：此处的α \alphaα不能是0或1。
Thm. 线性规划的多面体模型是凸集。
Def. 凸集的顶点：顶点无法表示成集合中其他元素的凸组合。
在这里插入图片描述
顶点的等价描述
从系数矩阵中抽取m列线性无关的列向量，组成可逆方阵。则由此可求得m个决策变量的值，再令其余的决策变量为0即可。
推论
顶点中正分量对应的系数向量线性无关。
一个线性规划问题标准模型最多有C n m C_{n}^{m}C
n
m
个顶点。
定义总结
基矩阵§：系数矩阵中抽取m列线性无关的列向量组成可逆方阵。
基本解：m个基变量有基矩阵和b ⃗ \vec{b}
b
决定，剩余(n-m)个变量都置0，称之为非基变量。
基本可行解（顶点）：基本解中可行的，即满足非负性约束
Thm. 线性规划标准模型的基本可行解就是可行集的顶点。
Thm. 标准模型的线性规划问题如有可行解，则定有基本可行解。
Thm. 线性规划标准模型中顶点的个数是有限的。
Thm. 线性规划标准模型的最优目标函数值如果有有限的目标函数值，则总在顶点处取到。
单纯形法
在顶点中沿着边搜索最优解的过程。
按照上述的原理，我们固然可以求出所有的基矩阵，进入求出所有的顶点。计算每一个顶点的目标函数值，找出其中最大的那个，但是这样做的计算量未免太大，因此有了单纯行法，即沿着边搜索顶点。
在这里插入图片描述
单纯形法就是一个不断的选择变量入基出基的过程。
假定已知一个基本可行解。（问题4）
如何计算选定进基变量后的基本可行解。(问题1）
如何选择进基变量使得目标函数值改善。（问题2）
如何判断已经找到最优的目标函数值。（问题3）
计算选定进基变量的基本可行解
Def. 基本可行解的表示式：基变量只出现在一个等式约束中。如：
在这里插入图片描述
此处的x 3 , x 4 , x 5 x_3,x_4,x_5x
3
,x
4
,x
5
就是基变量。
选定出基变量：保可行性的最小非负比值原理
由上所述，一个顶点对应一个基本可行解，其中m个基变量，(n-m)个非基变量。假定我们要选择某个非基变量x i x_ix
i
入基，实际上就是通过对增广矩阵做初等行变化使得x i x_ix
i
仅仅出现在一个等式约束中。比如我们通过变换，使得x i x_ix
i
仅仅出现在第j个等式约束中，如果此时仍然满足可行性，那么x i x_ix
i
就取代了原来在此处的基变量，成为新的基变量。
在进行初等行变换的过程中，要保证可行性，即
b ⃗ ≥ 0 \vec{b} \geq 0
b
≥0
。因此要选择最小非负比值。请看下面的例子：
在这里插入图片描述
假设我们要选择x 2 x_2x
2
入基，那么就是要通过初等行变换，使得x 2 x_2x
2
的系数向量中某一行是1，其余行都是0。如我们选择x 2 x_2x
2
仅出现在第3个等式约束中，即
在这里插入图片描述
则此时无法保证可行性，因为b ⃗ \vec{b}
b
中第1个分量是负数。
为了避免等式右侧出现负数，只能选择比值最小的一行，即第1行。即化成如下的形式：
在这里插入图片描述
如果此时我们想让x 3 x_3x
3
入基，此时的最小比值是第2行，即让该行为1，其余行为0。但是，为了让x 3 x_3x
3
的第二行为1，该行两端必须同时乘以一个负数，此时仍然无法保证b ⃗ ≥ 0 \vec{b} \geq0
b
≥0，因此只能选择系数非负的一行。
注：这里的非负性是指系数非负，而不是比值非负。即当b中某行分量是0，而该行入基变量系数是负数，仍不能入基。
在这里插入图片描述
特殊情况：没有非负比值，即没有有限的目标函数值。
在这里插入图片描述
选择入基变量的原则
选择某个入基变量使得目标函数能改善，通过检验数选择。
此处假设优化目标是求最大值。通过等式约束，将目标函数表示成非基变量的线性组合。即
f ( X ) = c 1 x j ( m + 1 ) + c 2 x j ( m + 2 ) + . . . + c n x j ( n ) + c o n s t f(X)=c_1x_{j(m+1)}+c_2x_{j(m+2)}+...+c_nx_{j(n)}+const
f(X)=c
1
x
j(m+1)
+c
2
x
j(m+2)
+...+c
n
x
j(n)
+const
只有选择检验数是正数的变量入基才有可能使得目标函数继续增大，因为入基之后变量只可能增大或者不变，而不可能减少。
如何确定已经找到了最优的目标函数值
此处假设优化目标是求最大值。
当每个非基变量的检验数都是负数时，目标函数已经达到了最大值。
退化情况
Thm. 收敛条件：每次迭代过程中，每个基本可行解的基变量都严格大于0，则每次迭代都能保证目标函数严格增加。而基本可行解的数目是有限的，因此上述过程不会一直进行下去，因此一定能在有限次迭代过程中找到最优解。
Def. 退化情况：某些基变量是0。则多个基矩阵对应同一个退化的顶点。
Thm. 循环迭代导致不收敛：多个基矩阵对应一个顶点，即每次出基入基都换了基矩阵，但是对应的退化顶点不变，即目标函数也不变。因此可能出现在几个基矩阵之间循环不止的情况。
避免退化：由于顶点的个数是有限的，我们只需标记那些已经迭代过的顶点，即可避免循环。
**bland法则：**始终选择下标最小的可入基和出基的变量。
当所有的基变量都严格大于0时，则这个基矩阵对应于非退化的顶点，此时可行基矩阵和顶点是一一对应的；
当某些基变量为0时，则这个基矩阵对应退化的顶点，一个退化的顶点对应数个可行基矩阵。
即给定一个可行基矩阵，一定能确定一个顶点，但是给定一个顶点时，其对应的基矩阵可能不唯一。
更一般地说，当顶点非退化时，可行基矩阵唯一；否则可行基矩阵不唯一。
如何确定初始的基本可行解
先将一般模型转化为标准模型，然后添加人工变量，在迭代过程中将人工变量都变成非基变量，则基变量就只剩下原来的变量。
在这里插入图片描述
大M法在这里插入图片描述
两阶段法
在这里插入图片描述
例题
本质就是不断的迭代单纯型表
在这里插入图片描述
在这里插入图片描述
一般线性规划问题总结
一般模型转化为标准型
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
基于单纯型表迭代的实质
求出非基变量的检验数
σ j ( k ) = c j ( k ) − C B T B − 1 P j ( k ) m + 1 ≤ k ≤ n \sigma_{j(k)}=c_{j(k)}-C_{B}^{T}B^{-1}P_{j(k)}\ m+1 \leq k \leq n
σ
j(k)
=c
j(k)
−C
B
T
B
−1
P
j(k)
m+1≤k≤n
确定进基变量
σ j ( t ) = m a x { σ j ( m + 1 ) , σ j ( m + 2 ) , . . . σ j ( n ) } \sigma_{j(t)} = max\{\sigma_{j(m+1)},\sigma_{j(m+2)},...\sigma_{j(n)}\}
σ
j(t)
=max{σ
j(m+1)
,σ
j(m+2)
,...σ
j(n)
}
确定出基变量
在这里插入图片描述
得到新的可行基矩阵
在这里插入图片描述
基于逆矩阵的单纯形法
在这里插入图片描述
核心问题：如何基于B − 1 B^{-1}B
−1
计算出B − 1 ~ \tilde{B^{-1}}
B
−1
~
。这两个矩阵仅仅有1列不一样，这是一个线性代数问题，与本课程的主要内容无关，此处不再赘述。
总结：单纯形法中可能遇到的3中特殊情况。
1. 退化问题：某些基变量为0
退化问题的现象是某些基变量为0，本质是一个退化的顶点对应多个可行基矩阵，后果是可能使得单纯形法不收敛。
在选择入基变量时，应该遵循blend法则，即每次选择可入基变量下标最小的那个。
2. 没有最小非负比值。
当选定入基变量后，需要根据“保证可行性的最小非负比值原理”选定出基变量，如果没有非负比值，则说明该变量可以趋于无穷，则该问题没有有限的最优目标函数值。
3. 某个非基变量的检验数为0.
在选择入基变量时，需要将目标函数表示成非基变量的表达式。以目标值是求最大问题的为例，此时应该选择检验数大于0的非基变量入基才能改善目标函数值。
当所有非基变量的检验数都为小于等于0的时候，无论选择谁入基，都会值得目标函数变得更差，因此这时候就达到了最优条件。
有一种特殊情况是某个非基变量的检验数为0，如果选取该变量入基，则目标函数值和原来一样，但是我们得到另一组不同的基本可行解，即最优目标函数值对应了多个基本可行解，这说明原问题有无穷多最优解。
4. 退化问题和非基变量检验数为0.
前者是一个顶点对应多个可行基矩阵，后者是最优目标函数值对应多个顶点。
前者可能导致单纯形法不收敛，后者说明该问题有无穷多解。
文章知识点与官方知识档案匹配
算法技能树首页概览
31789 人正在系统学习中
打开CSDN，阅读体验更佳
最优化技术——线性规划
最优化技术——线性规划 线性规划基本概念 线性规划问题就是在一组线性约束条件下，求解目标函数最优解的问题 标准形式 线性规划问题的标准形式： 目标函数求最大值 所有约束条件均由等式表示 每个约束条件右端常数常为非负值 所有决策变量为非负值 改造方法 所有的情况与改造方法 目标函数求最小值则应该改为求最大值： 方法——添加负号： minF=Σcjxj→maxF=−Σcjxj min F ...
继续访问
线性规划和对偶规划学习总结
在学习列生成和分解算法的时候，会频繁用到线性规划和对偶的知识，可以说LP和DP是基础。因此理解线性规划和对偶的本质是很重要的。 单纯形法是纯代数运算，从一个顶点跳到另一个顶点；且我们知道最优解一定在顶点取得，可是为什么在顶点一定会取得最优解？为什么从一个顶点跳到另一顶点，通过进基出基就可以实现，它俩为什么是一一对应的？除此以外，在分解算法中的极点和极射线和LP的解是什么样的对应关系？ 这些问题，应该说必须了解了线性规划和对偶的本质才能够深刻理解，而不至于陷入机械无脑的计算。大概看了慕课的课程，感觉山东大
继续访问
最新发布 03 线性规划模型
03 线性规划模型
继续访问
第五章 线性规划方法 Linear Programming
第五章 线性规划方法 Linear Programming5.1 线性规划问题的一般形式5.2 线性规划问题的解5.2.1 基本解的产生与转换5.2.2 基本可行解的产生与转换5.2.3 基本可行解的变换条件1. 最优性条件2. 非负性条件5.3 单纯形算法 The Simplex Method 5.1 线性规划问题的一般形式 5.2 线性规划问题的解 基本解： 只满足约束方程的解。 基本可行解： 同时满足约束方程和变量非负约束的解。 最优解： 使目标函数取得最小值的基本可行解。 5.2.1 基本解的产生与
继续访问
关于数学建模中线性规划总结
一、python方法解决 from scipy import optimize as op import numpy as np c=np.array([2,3,-5]) c = np.array([2,3,-5]) A = np.array([[-2,5,-1],[1,3,1]]) b= np.array([-10,12]) Aeq = np.array([[1,1,1]]) beq = np.array([7]) #求解 res = op.linprog(-c,A,b,Aeq,beq) print(
继续访问
八、线性规划 顶点、极值点和基本可行解决方案
假设我们正在求解方程形式的一般线性程序： 这里，是一个的矩阵，，，今天，我们将假设 的行是线性独立的。 （如果不是，那么系统 没有解，或者某些方程是多余的。在第一种情况下，我们只是忘记分析这样的线性程序；在第二种情况下，我们可以从删除冗余行。） 我们已经非正式地说过，基本可行的解决方案是“尽可能多的变量”为0。这不是很精确：在某些情况下（由于退化），可能有异常多的0值，并且我们不希望这与我们的定义混淆。 相反，我们进行如下定义。 选择一些列（或变量） 的 做为
继续访问
【算法设计zxd】第3章迭代法04 线性规划
线性规划 研究线性约束条件下线性目标函数 的极值问题的数学理论和方法。 线性规划问题形式化表达 目标函数 约束条件 线性规划问题的可行性解 线性规划问题的可行区域 线性规划问题的最优解（x1,x2,……，xn的值） 线性规划问题的最优值 单纯形算法特点 (1) 只对约束条件的若干组合进行测试，测试的毎一步都使 目标函数的值向期望值逼近； (2) 一般经过不大于m或n次迭代就可求得最优解。 线性规划标准形式 (1)它必须是一个最大化问题。如果是..
继续访问
线性规划部分概念及重要性质（运筹学导论笔记）
模型解的术语 可行解：满足所有约束条件的解 非可行解：至少一个约束条件不被满足的解 可行域：所有可行解的集合 最优解：目标函数取得最有价值的可行解 顶点可行解（CPF）：位于可行域顶点的解 顶点可行解与最优解的关系：考虑任意具有可行解与有界可行域的线性规划问题，一定具有顶点可行解和至少一个最优解，而且，最优的顶点可行解一定是最优解；因此，若一个问题恰有一个最优解，它一定是顶点可行解，若一个问题有多个最优解，其中至少两个一定是顶点可行解 比例性假设：每个活动对于目标函数值Z的贡献与活动级别xj成比例的 可加性
继续访问
Mathematics for Machine Learning--学习笔记(线性无关)
1.5 Linear Independence（线性无关） 接下来就要学习如何处理向量了。首先，我们先介绍线性组合和线性无关的概念。 Linear Combination(线性组合)：存在一个向量空间V和有限的x1,⋯ ,xk∈Vx_1,\cdots,x_k\in Vx1,⋯,xk∈V.每一个元素vvv都有如下形式：v=λ1x1+⋯+λkxk=∑i=1kλixi∈Vv=\lambda_1 x_1+\cdots+\lambda_k x_k=\sum_{i = 1}^{k} {\lambda_i x_i
继续访问
线性规划——规范型，标准型，基阵、基本解、基本可行解、基变量、非基变量.... 概念梳理
文章目录前言最优化—线性规划模型问题线性规划模型的一般形式（min）线性规划规范形式线性规划标准型模型的转换线性规划中的规律规范形式顶点的数学描述标准形式顶点的数学描述标准形式顶点的等价描述之一标准形式顶点的等价描述之二线性规划标准形式的一些基本概念线性规划标准形式的基本定理 前言 此总结参考 清华 王焕刚老师的教程。 最优化—线性规划 模型问题 线性规划模型的一般形式（min） min∑j=1ncjxj s.t. ∑j=1naijxj=bi,∀1≤i≤p∑j=1naijxj≥bi,∀
继续访问
最优化——线性规划总结1(线性规划标准型，规范型，顶点)
线性规划的形式 标准型 规范型 线性规划的求解思路 前提条件 线性规划：凸优化（凸集上的凸函数的优化） 线性规划的可行集是凸集，优化函数是凸函数（仿射函数嘛） 总有顶点是最优解，所有顶点组成的集合总是有限集，所以可以在顶点集中找到最优解。 主要思路 根据前提条件来看，我们求解线性规划的思路：找到所有的顶点，在顶点中找到最优的那个，就是最优解。相当于缩小了搜索范围。 怎么搞 首先计算顶点：顶点是改点所有起作用约束构成的线性方程组的唯一解。 因为所有的线性规划形式都能转换成标准型，所以这里只考虑标准型的
继续访问
线性规划图解法求最优解_高考数学【线性规划】知识点相关解析~
一、知识梳理1、目标函数：P=2x＋y是一个含有两个变量x和y的函数，称为目标函数。2、可行域：约束条件表示的平面区域称为可行域。3、整点：坐标为整数的点叫做整点。4、线性规划问题：求线性目标函数在线性约束条件下的最大值或最小值的问题，通常称为线性规划问题。只含有两个变量的简单线性规划问题可用图解法来解决。5、整数线性规划：要求量整数的线性规划称为整数线性规划。二、疑难知识导析线性规划是...
继续访问
算法最优化（2）线性规划问题中的常见概念辨析：可行解，最优解，基，基向量，非基向量，基变量，非基变量等等
zz
继续访问
【线性规划】基本概念
线性规划的概念 线性规划(Linear Programming 简记 LP)是了运筹学中数学规划的一个重要分支。自从 1947 年 G. B. Dantzig 提出 求解线性规划的单纯形法以来，线性规划在理论上趋向成熟，在实用中由于计算机能处理成千上万个约束条件和决策变量的线性规划问题之后，线性规划现代管理中经常采用的基本方法之一。 在解决实际问题时，需要把问题归结成一个线性规划数学模型，关键及难点在于选适当的决策变量建立恰当的模型，这直接影响到问题的求解。 线性规划问题的目标函数及约束条件均为线性函数；约
继续访问
【运筹学】什么是基变量？对于线性规划问题中“基”概念的理解（3月3日学习笔记）
在学习《线性规划与目标规划》的过程中，课程的主讲老师郭韧给出了对于基概念的定义如下图 图片来源:运筹学（中国大学mooc网） 由此我产生了几个疑惑:1.如何理解B是线性规划问题的一个基？2.为什么说最多有CnmC_n^mCnm个基呢？ 1.如何理解B是线性规划问题的一个基？1.如何理解B是线性规划问题的一个基？1.如何理解B是线性规划问题的一个基？ 在回答第一个...
继续访问
【运筹学】线性规划 最优解分析 ( 唯一最优解 | 无穷多最优解 | 无界解 | 无可行解 | 迭代范围 | 求解步骤 )
一、唯一最优解、 二、无穷多最优解、 三、无界解、 四、无可行解、 五、线性规划迭代范围、 六、线性规划求解步骤
继续访问
线性规划与非线性规划的求解
单纯形法求解线性规划 一、大M法求解线性规划的原理 （1）、大M法首先将线性规划问题化为标准型。如果约束方程组中包含有一个单位矩阵I，那么已经得到了一个初始可行基。否则在约束方程组的左边加上若千个非负的人工变量，使人工变量对应的系数列向量与其它变量的系数列向量共同构成-一个单位矩阵。以单位矩阵为初始基，即可求得一-个初始的基本可行解。 为了求得原问题的初始基本可行解，必须尽快通过迭代过程把人工变量...
继续访问
热门推荐 线性规划算法详解
线性规划 首先什么是线性规划，大致的定义我总结为在线性的目标和约束中，找出一个最优解。 举个例子： M1和M2两种原料用于生产内外墙涂料，M1日最大可用量24吨，M2日最大可用量为6吨，外墙涂料每吨需要6吨M1，1吨M2，内墙涂料每吨需要4吨M12，吨M2，外墙涂料每吨利润5个单位，内墙涂料每吨利润4个单位。且市场需求调查数据得出，内墙日需求量不超过外墙的日需求量+1吨，内墙最大日需求量为...
继续访问
运筹学 —线性规划总结
线性规划问题 1. 概述 线性规划问题是在一组线性约束下，求资源配置的最大最小值的问题。 直观的变现是在一个约束条件围成的区域上寻找一个点，这个点使得资源配置最优化： 2. 线性规划的思想 线性规划的思路是将不等式转换为等式，最终求得一个满足等式的解。 下面的约束式必然可以转换为[P|N]*X=B的形式，这里P是线性无关的M*M的方正。
继续访问
最优化——退化和某个非基变量检验数为零
文章目录退化和某个非基变量检验数为零退化问题退化问题的本质某个非基变量检验数为零 退化和某个非基变量检验数为零 退化问题 基本可行解的基变量数值等于0。 退化问题的本质 多个可行基阵对应于一个基本可行解。 某个非基变量检验数为零 对于求max的线性规划问题，如果所有检验数均满足 则说明已经得到最优解， 若此时某非基变量的检验数为零 ，则说明该优化问题有无穷多最优解。 ...
如何根据概率密度函数生成随机数
我这里并不是要讲“伪随机”、“真随机”这样的问题，而是关于如何生成服从某个概率分布的随机数（或者说 sample）的问题。比如，你想要从一个服从正态分布的随机变量得到 100 个样本，那么肯定抽到接近其均值的样本的概率要大许多，从而导致抽到的样本很多是集中在那附近的。当然，要解决这个问题，我们通常都假设我们已经有了一个 生成 0 到 1 之间均匀分布的随机数的工具，就好像 random.org 给我们的结果那样，事实上许多时候我们也并不太关心它们是真随机数还是伪随机数，看起来差不多就行了。 :p
现在再回到我们的问题，看起来似乎是很简单的，按照概率分布的话，只要在概率密度大的地方多抽一些样本不就行了吗？可是具体要怎么做呢？要真动起手 来，似乎有不是那么直观了。实际上，这个问题曾经也是困扰了我很久，最近又被人问起，那我们不妨在这里一起来总结一下。为了避免一下子就陷入抽象的公式推 导，那就还是从一个简单的具体例子出发好了，假设我们要抽样的概率分布其概率密度函数为 p(x) = \frac{1}{9}x^2 ，并且被限制在区间 [0, 3] 上，如右上图所示。
好了，假设现在我们要抽 100 个服从这个分布的随机数，直观上来讲，抽出来的接近 3 的数字肯定要比接近 0 的数字要多。那究竟要怎样抽才能得到这样的结果呢？由于我们实际上是不能控制最原始的随机数生成过程的，我们只能得到一组均匀分布的随机数，而这组随机数 的生成过程对于我们完全是透明的，所以，我们能做的只有把这组均匀分布的随机数做一些变换让他符合我们的需求。找到下手的点了，可是究竟要怎样变换呢？有 一个变换相信大家都是很熟悉的，假设我们有一组 [0,1] 之间的均匀分布的随机数 X_0 ，那么令 X_1=3X_0 的话，X_1 就是一组在 [0,3] 之间均匀分布的随机数了，不难想象，X_1 等于某个数 x^* 的概率就是 X_0 等于 x^*/3 的概率（“等于某个数的概率”这种说法对于连续型随机变量来说其实是不合适的，不过大概可以理解所表达的意思啦）。似乎有一种可以“逆转回去”的感觉了。
于是让我们来考虑更一般的变换。首先，我们知道 X_1 的概率密度函数是 f(x) = 1/3, x\in[0,3] ，假设现在我们令 Y = \phi (X_1) ，不妨先假定 \phi(\cdot) 是严格单调递增的函数，这样我们可以求其逆函数 \phi^{-1}(\cdot) （也是严格单调递增的）。现在来看变换后的随机变量 Y 会服从一个什么样的分布呢？
这里需要小心，因为这里都是连续型的随机变量，并不像离散型随机变量那样可以说成“等于某个值的概率”，因此我们需要转换为概率分布函数来处理，也就是求一个积分啦：
\displaystyle F(x) = P(X \leq x) = \int_{-\infty}^x f(t)dt
那么 X_1 的概率分布函数为 F(x) = \frac{1}{3}x 。很显然 Y 小于或等于某个特定的值 y^* 这件事情是等价于 X_1=\phi^{-1}(Y)\leq\phi^{-1}(y^*) 这件事情的。换句话说，P(Y\leq y^*) 等于 P(X_1 \leq \phi^{-1}(y^*)) 。于是，Y 的概率分布函数就可以得到了：
\displaystyle G(y) = P(Y \leq y) = P(X_1 \leq \phi^{-1}(y)) = F(\phi^{-1}(y))
再求导我们就能得到 Y 的概率密度函数：
\displaystyle g(y) = \frac{dG(y)}{dy} = f(\phi^{-1}(y))\frac{d}{dy}\phi^{-1}(y)
这样一来，我们就得到了对于一个随机变量进行一个映射 \phi(\cdot) 之后得到的随即变量的分布，那么，回到我们刚才的问题，我们想让这个结果分布就是我们所求的，然后再反推得 \phi(\cdot) 即可：
\displaystyle \frac{1}{9}y^2 = g(y) = f(\phi^{-1}(y))\frac{d}{dy}\phi^{-1}(y) = \frac{1}{3}\frac{d}{dy}\phi^{-1}(y)
经过简单的化简就可以得到 \phi^{-1}(y) = \frac{1}{9} y^3 ，亦即 \phi(x) = (9x)^{1/3} 。也就是说，把得到的随机数 X_1 带入到到函数 \phi(\cdot) 中所得到的结果，就是符合我们预期要求的随机数啦！ :D 让我们来验证一下：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
#!/usr/bin/python import numpy as np import matplotlib.pyplot as plot N = 10000 X0 = np.random.rand(N) X1 = 3*X0 Y = np.power(9*X1, 1.0/3) t = np.arange(0.0, 3.0, 0.01) y = t*t/9 plot.plot(t, y, 'r-', linewidth=1) plot.hist(Y, bins=50, normed=1, facecolor='green', alpha=0.75)plot.show()
这就没错啦，目的达成啦！让我们来总结一下。问题是这样的，我们有一个服从均匀分布的随机变量 X ，它的概率密度函数为一个常数 f(x)=C ，如果是 [0,1] 上的分布，那么常数 C 就直接等于 1 了。现在我们要得到一个随机变量 Y 使其概率密度函数为 g(y) ，做法就是构造出一个函数 \phi(\cdot) 满足（在这里加上了绝对值符号，这是因为 \phi(\cdot) 如果不是递增而是递减的话，推导的过程中有一处就需要反过来）
\displaystyle g(y) = f(\phi^{-1}(y))\left|\frac{d}{dy}\phi^{-1}(y)\right| = C\left|\frac{d}{dy}\phi^{-1}(y)\right|
反推过来就是，对目标 y 的概率密度函数求一个积分（其实就是得到它的概率分布函数 CDF ，如果一开始就拿到的是 CDF 当然更好），然后求其反函数就可以得到需要的变换 \phi(\cdot) 了。实际上，这种方法有一个听起来稍微专业一点的名字：Inverse Transform Sampling Method 。不过，虽然看起来很简单，但是实际操作起来却比较困难，因为对于许多函数来说，求逆是比较困难的，求积分就更困难了，如果写不出解析解，不得已只能用数 值方法来逼近的话，计算效率就很让人担心了。可事实上也是如此，就连我们最常见的一维标准正态分布，也很难用这样的方法来抽样，因为它的概率密度函数
\displaystyle g(y) = \frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}y^2}
的不定积分没有一个解析形式。这可真是一点也不好玩，费了这么大劲，结果好像什么都干不了。看来这个看似简单的问题似乎还是比较复杂的，不过也不要灰心，至少对于高斯分布来说，我们还有一个叫做 Box Muller 的方法可以专门来做这个事情。因为高斯分布比较奇怪，虽然一维的时候概率分布函数无法写出解析式，但是二维的情况却可以通过一些技巧得出一个解析式来。
首先我们来考虑一个二维的且两个维度相互独立的高斯分布，它的概率密度函数为
\displaystyle f(x,y) = \frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}\cdot\frac{1}{\sqrt{2\pi}}e^{-\frac{y^2}{2}} = \frac{1}{2\pi}e^{-\frac{x^2+y^2}{2}}
这个分布是关于原点对称的，如果考虑使用极坐标 (\theta,r) （其中 \theta\in[0,2\pi), r\in[0,\infty) ）的话，我们有 x = r\cos\theta，y=r\sin\theta 这样的变换。这样，概率密度函数是写成：
\displaystyle f(\theta,r) = \frac{1}{2\pi}e^{-\frac{r^2}{2}}
注意到在给定 r 的情况下其概率密度是不依赖于 \theta 的，也就是说对于 \theta 来说是一个均匀分布，这和我们所了解的标准正态分布也是符合的：在一个圆上的点的概率是相等的。确定了 \theta 的分布，让我们再来看 r，用类似于前面的方法：
\displaystyle \begin{aligned} P(rR) = \int_0^{2\pi}\int_0^R\frac{1}{2\pi}e^{\frac{r^2}{2}}rdrd\theta \ = \int_0^Re^{-\frac{r^2}{2}}rdr \ = 1-e^{-\frac{R^2}{2}} \end{aligned}
根据前面得出的结论，我现在得到了 r 的概率分布函数，是不是只要求一下逆就可以得到一个 \phi(\cdot) 了？亦即 \phi(t) = \sqrt{-2\log (1-t)} 。
现在只要把这一些线索串起来，假设我们有两个相互独立的平均分布在 [0,1] 上的随机变量 T_1 和 T_2 ，那么 2\pi T_1 就可以得到 \theta 了，而 \phi(T_2) = \sqrt{-2\log(1-T_2)} 就得到 r 了（实际上，由于 T_2 和 1-T_2 实际上是相同的分布，所以通常直接写为 \sqrt{-2\log T_2}）。再把极坐标换回笛卡尔坐标：
\displaystyle \begin{aligned} x = r\cos\theta = \sqrt{-2\log T_2}\cdot \cos(2\pi T_1) \ y = r\sin\theta = \sqrt{-2\log T_2}\cdot \sin(2\pi T_1) \end{aligned}
这样我们就能得到一个二维的正态分布的抽样了。可以直观地验证一下，二维不太好画，就画成 heatmap 了，看着比较热的区域就是概率比较大的，程序如下：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
#!/usr/bin/python import numpy as np import matplotlib.pyplot as plot N = 50000 T1 = np.random.rand(N) T2 = np.random.rand(N) r = np.sqrt(-2*np.log(T2)) theta = 2*np.pi*T1 X = r*np.cos(theta) Y = r*np.sin(theta) heatmap, xedges, yedges = np.histogram2d(X, Y, bins=80) extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]] plot.imshow(heatmap, extent=extent) plot.show()
画出来的图像这个样子：
不太好看，但是大概的形状是可以看出来的。其实有了二维的高斯分布，再注意到两个维度在我们这里是相互独立的，那么直接取其中任意一个维度，就是一个一维高斯分布了。如下：
如果 X\sim N(0,1) 即服从标准正态分布的话，则有 \sigma X+\mu \sim N(\mu, \sigma^2) ，也就是说，有了标准正态分布，其他所有的正态分布的抽样也都可以完成了。这下总算有点心满意足了。不过别急，还有最后一个问题：多元高斯分布。一般最常 用不就是二元吗？二元不是我们一开始就推出来了吗？推出来了确实没错，不过我们考虑的是最简单的情形，当然同样可以通过 \sigma X+\mu 这样的方式来处理每一个维度，不过高维的情形还有一个需要考虑的就是各个维度之间的相关性——我们之前处理的都是两个维度相互独立的情况。对于一般的多维正态分布 X\sim N(\mathbf{\mu}, \Sigma) ，如果各个维度之间是相互独立的，就对应于协方差矩阵 \Sigma 是一个对角阵，但是如果 \Sigma 在非对角线的地方存在非零元素的话，就说明对应的两个维度之间存在相关性。
这个问题还是比较好解决的，高斯分布有这样的性质：类似于一维的情况，对于多维正态分布 X\sim N(\mathbf{\mu}, \Sigma)，那么新的随机变量 X_1=\mathbf{\mu}_1 + LX 将会满足
\displaystyle X_1 \sim N(\mathbf{\mu}_1+L\mu, L\Sigma L^T)
所以，对于一个给定的高斯分布 N(\mathbf{\mu}, \Sigma) 来说，只要先生成一个对应维度的标准正态分布 X\sim N(0, I) ，然后令 X_1 = \mu+LX 即可，其中 L 是对 \Sigma 进行 Cholesky Decomposition 的结果，即 \Sigma = LL^T 。
结束之前让我们来看看 matlab 画个 3D 图来改善一下心情：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
N = 50000; T1 = rand(1, N); T2 = rand(1, N); r = sqrt(-2*log(T2)); theta = 2*pi*T1; X =[r.*cos(theta); r.*sin(theta)]; mu = [1; 2]; Sigma = [5 2; 2 1]; L = chol(Sigma); X1 = repmat(mu,1, N) + L*X; nbin = 30; hist3(X1', [nbin nbin]); set(gcf, 'renderer', 'opengl'); set(get(gca,'child'), 'FaceColor', 'interp', 'CDataMode', 'auto'); [z c] = hist3(X1', [nbin nbin]); [x y] =meshgrid(c{1}, c{2}); figure; surfc(x,y,-z);
下面两幅图，哪幅好看一些（注意坐标比例不一样，所以看不出形状和旋转了）？似乎都不太好看，不过感觉还是比前面的 heatmap 要好一点啦！
然后，到这里为止，我们算是把高斯分布弄清楚了，不过这只是给一个介绍性的东西，里面的数学推导也并不严格，而 Box Muller 也并不是最高效的高斯采样的算法，不过，就算我们不打算再深入讨论高斯采样，采样这个问题本身也还有许多不尽人意的地方，我们推导出来的结论可以说只能用 于一小部分简单的分布，连高斯分布都要通过 trick 来解决，另一些本身连概率密度函数都写不出来或者有各种奇怪数学特性的分布就更难处理了。所以本文的标题里也说了，这是上篇，如果什么时候有机会抽出时间 来写下篇的话，我将会介绍一些更加通用和强大的方法，诸如 Rejection Sampling 、Gibbs Sampling 以及 Markov Chain Monte Carlo (MCMC) 等方法。如果你比较感兴趣，可以先自行 Google 一下解馋！ :D